module atmos_tracer_driver_mod
! <CONTACT EMAIL="William.Cooke@noaa.gov">
!   William Cooke
! </CONTACT>

! <REVIEWER EMAIL="Matthew.Harrison@noaa.gov">
!   Matt Harrison
! </REVIEWER>

! <REVIEWER EMAIL="Bruce.Wyman@noaa.gov">
!   Bruce Wyman
! </REVIEWER>

! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/"/>

! <OVERVIEW>
!     This code allows the user to easily add tracers to the FMS framework.
! </OVERVIEW>

! <DESCRIPTION>
!
!     This code allows a user to easily implement tracer code in the FMS
!     framework.  The tracer and tracer tendency arrays are supplied along
!     with longtitude,  latitude, wind, temperature, and pressure
!     information which allows a  user to implement sources and sinks of the
!     tracer which depend on these parameters.
!
!     In the following example, radon being implemented in the atmosphere
!     will be used as an example of how to implement a tracer in the FMS
!     framework.
!
!     Within the global scope of tracer_driver_mod a use
!     statement should be inserted for each tracer to be added.
!<PRE>      use radon_mod, only : radon_sourcesink, radon_init, radon_end </PRE>
!
!     An integer parameter, which  will be used as an identifier for the
!     tracer, should be assigned.
!<PRE>
!      integer :: nradon
!</PRE>
!     Within tracer_driver_init a call to the tracer manager is needed in
!     order  to identify which tracer it has set the tracer as.
!<PRE>
!      nradon = get_tracer_index(MODEL_ATMOS,'radon')
!</PRE>
!     Here MODEL_ATMOS is a parameter defined in field_manager.
!          'radon' is the name of the tracer within the field_table.
!
!     If the tracer exists then the integer returned will be positive and it
!     can be used to call the initialization routines for the individual
!     tracers.
!<PRE>
!      if (nradon > 0) then
!           call radon_init(Argument list)
!      endif
!</PRE>
!
!     Within tracer_driver the user can also use the identifier to surround
!     calls to the source-sink routines for the tracer of interest.
!
!<PRE>
!      if (nradon > 0 .and. nradon <= nt) then
!          call radon_sourcesink (Argument list)
!          rdt(:,:,:,nradon)=rdt(:,:,:,nradon)+rtnd(:,:,:)
!      endif
!</PRE>
!
!     It is the users responsibility to add the tendency generated by the
!     sourcesink routine.
!
!     Within tracer_driver_end the user can add calls to the
!     terminators for the appropriate source sink routines.
!
!<PRE>      call radon_end</PRE>
!
!     This may simply be a deallocation statement or a routine to send
!     output to the logfile stating that the termination routine has been
!     called.
!
! <NOTE>
! This code has been modified by Paul.Ginoux@noaa.gov in 2005
! to include the following aerosols:
!  - SO4 with a simplified SOx chemistry
!  - sea salt with a size distributon splitted into 5 bins
!  - mineral dust with a size distributon splitted into 5 bins
!  - carbonaceous aerosols with new emissions by Shekar.Reddy@noaa.gov
! </NOTE>
!
! </DESCRIPTION>


!-----------------------------------------------------------------------

use mpp_mod,               only : input_nml_file
use fms_mod,               only : file_exist, close_file,&
                                  open_namelist_file, check_nml_error, &
                                  write_version_number, &
                                  error_mesg, &
                                  FATAL, &
                                  mpp_pe, &
                                  mpp_root_pe, &
                                  stdlog, stdout, &
                                  mpp_clock_id, &
                                  mpp_clock_begin, &
                                  mpp_clock_end, &
                                  CLOCK_MODULE, &
                                  uppercase
use time_manager_mod,      only : time_type, &
                                  get_date, get_date_julian, &
                                  real_to_time_type
use diag_manager_mod,      only : register_diag_field, send_data
use atmos_cmip_diag_mod,   only : register_cmip_diag_field_3d, &
                                  register_cmip_diag_field_2d, &
                                  send_cmip_data_3d, &
                                  cmip_diag_id_type, &
                                  query_cmip_diag_id

use astronomy_mod,         only : astronomy_init, diurnal_solar, universal_time
use tracer_manager_mod,    only : get_tracer_index,   &
                                  get_number_tracers, &
                                  get_tracer_names,   &
                                  get_tracer_indices, &
                                  adjust_positive_def, &
                                  query_method, &
				  NO_TRACER
use field_manager_mod,     only : MODEL_ATMOS
use atmos_tracer_utilities_mod, only :                      &
                                  dry_deposition,           &
                                  dry_deposition_init,      &
                                  dry_deposition_time_vary, &
                                  dry_deposition_endts,     &
                                  atmos_tracer_utilities_init, &
                                  get_rh, get_w10m, get_cldf, &
                                  sjl_fillz, &
                                  get_cmip_param, get_chem_param
use constants_mod,         only : grav, WTMAIR, PI, AVOGNO
use atmos_radon_mod,       only : atmos_radon_sourcesink,   &
                                  atmos_radon_init,         &
                                  atmos_radon_end
use atmos_carbon_aerosol_mod, only : &
                                  atmos_carbon_aerosol_time_vary,  &
                                  atmos_carbon_aerosol_endts,      &
                                  atmos_carbon_aerosol_driver,  &
                                  atmos_carbon_aerosol_init,&
                                  atmos_carbon_aerosol_end
use atmos_convection_tracer_mod,only: &
                                  atmos_convection_tracer_init, &
                                  atmos_cnvct_tracer_sourcesink, &
                                  atmos_convection_tracer_end
use atmos_sulfur_hex_mod,  only : atmos_sf6_sourcesink,     &
                                  atmos_sulfur_hex_init,    &
                                  atmos_sulfur_hex_end
use atmos_ch3i_mod,        only : atmos_ch3i, &
                                  atmos_ch3i_init, &
                                  atmos_ch3i_time_vary, &
                                  atmos_ch3i_endts, &
                                  atmos_ch3i_end
use atmos_sea_salt_mod,    only : atmos_sea_salt_sourcesink,     &
                                  atmos_sea_salt_init,    &
                                  atmos_sea_salt_end, &
                                  is_seasalt_tracer,          &
                                  do_seasalt
use atmos_dust_mod,        only : atmos_dust_sourcesink,   &
                                  atmos_dust_init,         &
                                  atmos_dust_flux_init,    &
                                  atmos_dust_gather_data,  &
                                  atmos_dust_time_vary,    &
                                  atmos_dust_endts,        &
                                  is_dust_tracer,          &
                                  dust_has_surf_setl_flux, &
                                  get_dust_surf_setl_flux, &
                                  atmos_dust_end,          &
                                  do_dust
use atmos_sulfate_mod,     only : atmos_sulfate_init, &
                                  atmos_sulfate_time_vary, &
                                  atmos_sulfate_endts,     &
                                  atmos_sulfate_end, &
                                  atmos_DMS_emission, &
                                  atmos_SOx_emission, &
                                  atmos_SOx_chem
use atmos_soa_mod,         only : atmos_SOA_init, &
                                  atmos_SOA_time_vary, &
                                  atmos_SOA_endts, &
                                  atmos_SOA_end, &
                                  atmos_SOA_chem
use tropchem_driver_mod,   only : tropchem_driver, &
                                  tropchem_driver_time_vary, &
                                  tropchem_driver_endts, &
                                  tropchem_driver_init
use atmos_regional_tracer_driver_mod,only : regional_tracer_driver, &
                                            regional_tracer_driver_init
use strat_chem_driver_mod, only : strat_chem, strat_chem_driver_init
use atmos_age_tracer_mod,  only : atmos_age_tracer_init, atmos_age_tracer, &
                                  atmos_age_tracer_end
use atmos_co2_mod,         only : atmos_co2_sourcesink,   &
                                  atmos_co2_emissions,          &
                                  atmos_co2_time_vary,          &
                                  atmos_co2_gather_data,        &
                                  atmos_co2_flux_init,          &
                                  atmos_co2_init,               &
                                  atmos_co2_end
use atmos_tropopause_mod,only: &
                                  atmos_tropopause_init, &
                                  atmos_tropopause

use interpolator_mod,      only : interpolate_type

implicit none
private
!-----------------------------------------------------------------------
!----- interfaces -------

public  atmos_tracer_driver,            &
        atmos_tracer_driver_init,       &
        atmos_tracer_driver_time_vary,       &
        atmos_tracer_driver_endts,       &
        atmos_tracer_driver_end,        &
        atmos_tracer_flux_init,         &
        atmos_tracer_driver_gather_data, &
        atmos_tracer_driver_gather_data_down, &
        atmos_tracer_has_surf_setl_flux, &
        get_atmos_tracer_surf_setl_flux

!-----------------------------------------------------------------------
!----------- namelist -------------------
logical :: prevent_flux_through_ice = .false.  , step_update_tracer = .false.
                               ! when true, tracers will only be fluxed
                               ! through the non-ice-covered portions of
                               ! ocean grid boxes


namelist /atmos_tracer_driver_nml / prevent_flux_through_ice, step_update_tracer
!-----------------------------------------------------------------------
!
!  When initializing additional tracers, the user needs to make the
!  following changes.
!
!  Add an integer variable below for each additional tracer.
!  This should be initialized to zero.
!
!-----------------------------------------------------------------------

! Indices for timing the various chemistry routines
integer :: radon_clock = 0
integer :: convect_clock = 0
integer :: age_tracer_clock = 0
integer :: stratozone_clock = 0
integer :: tropchem_clock = 0
integer :: carbon_clock = 0
integer :: dust_clock = 0
integer :: seasalt_clock = 0
integer :: sulfur_clock = 0
integer :: SOA_clock = 0
integer :: sf6_clock = 0
integer :: ch3i_clock = 0
integer :: co2_clock = 0
integer :: tropopause_clock = 0

logical :: do_tropchem = .false.  ! Do tropospheric chemistry?
logical :: do_coupled_stratozone = .FALSE. !Do stratospheric chemistry?

integer :: nsphum  ! Specific humidity parameter

integer :: no3 = 0
integer :: no3ch = 0
integer :: nextinct = 0
integer :: naerosol = 0
integer :: nbcphobic =0
integer :: nbcphilic =0
integer :: nomphobic =0
integer :: nomphilic =0
integer :: nclay     =0
integer :: nsilt     =0
integer :: nsf6      =0
integer :: nDMS      =0
integer :: nSO2      =0
integer :: nSO4      =0
integer :: nMSA      =0
integer :: nSOA      =0
integer :: nH2O2     =0
integer :: nch3i     =0
integer :: nage      =0
integer :: nco2      =0
integer :: nNH4NO3   =0
integer :: nNH4      =0
integer :: nDMS_cmip =0
integer :: nSO2_cmip =0
integer :: nSO4_cmip =0
integer :: nNH3_cmip =0
integer :: nOH       =0
integer :: nC4H10    =0
integer :: ncodirect =0
integer :: ne90 =0
integer :: nsulfate  =0
integer, dimension(5) :: tr_nbr_sulfate=0
logical, dimension(5) :: do_tracer_sulfate=.false.

real    :: ozon(11,48),cosp(14),cosphc(48),photo(132,14,11,48),   &
           solardata(1801),chlb(90,15),ozb(144,90,12),tropc(151,9),  &
           dfdage(90,48,8),anoy(90,48)


integer, dimension(:), pointer :: nradon
integer, dimension(:), pointer :: nconvect

integer :: nt     ! number of activated tracers
integer :: ntp    ! number of activated prognostic tracers

logical :: use_tau=.false.

character(len=6), parameter :: module_name = 'tracer'
character(len=7), parameter :: mod_name = 'tracers'

logical :: module_is_initialized = .FALSE.

type(interpolate_type), allocatable :: drydep_data(:)


integer, allocatable :: local_indices(:)
! This is the array of indices for the local model.
! local_indices(1) = 5 implies that the first local tracer is the fifth
! tracer in the tracer_manager.

!<f1p
integer, dimension(:), allocatable :: id_tracer_diag
integer, dimension(:,:), allocatable :: id_tracer_diag_hour
integer :: id_ps_hour(24),id_temp_hour(24),id_local_hour
!>
integer :: id_landfr, id_seaicefr, id_snowfr, id_vegnfr, id_vegnlai
integer :: id_om_ddep, id_bc_ddep, &
           id_nh4_ddep_cmip
integer :: id_nh4no3_col, id_nh4_col
integer :: id_nh4no3_cmip, id_nh4_cmip
integer :: id_nh4no3_cmipv2, id_nh4_cmipv2
integer :: id_so2_cmip, id_dms_cmip
integer :: id_so2_cmipv2, id_dms_cmipv2
integer :: id_n_ddep, id_n_ox_ddep, id_n_red_ddep

 type(cmip_diag_id_type) :: ID_concno3, ID_concnh4, ID_concso2, ID_concdms
 type(cmip_diag_id_type) :: ID_airmass, ID_pm1, ID_pm10, ID_pm25, ID_OM, ID_BC, ID_DUST, ID_SS

 integer :: id_sconcno3, id_sconcnh4, id_loadno3, id_loadnh4
 integer :: id_dryso2, id_dryso4, id_drydms, id_drynh3, &
            id_drynh4, id_drybc, id_drypoa, id_drysoa, id_dryoa

!for cmip6 (f1p)
 type(cmip_diag_id_type), allocatable :: ID_tracer_mol_mol(:),ID_tracer_kg_kg(:)
 integer, allocatable :: id_tracer_surf_mol_mol(:),id_tracer_surf_kg_kg(:)
 integer, allocatable :: id_tracer_col_kg_m2(:)
 integer              :: id_pm25_surf, id_dust_col_kg_m2, id_seasalt_col_kg_m2, id_bc_col_kg_m2
 integer              :: id_toz, id_tropoz
 real, allocatable    :: conv_vmr_mmr(:), nb_N(:),nb_N_ox(:),nb_N_red(:), frac_pm1(:), frac_pm10(:), frac_pm25(:)
 integer, allocatable :: id_tracer_ddep_kg_m2_s(:)

 real, parameter      :: o3_column_factor = 2.687e25 ! Molecular density (molec/m3) of air at STP (T=0C, P=1atm)
                                                     ! Used to convert from molec/m2 to equivalent depth (in m) at STP

!-----------------------------------------------------------------------
type(time_type) :: Time

!---- version number -----
character(len=128) :: version = '$Id$'
character(len=128) :: tagname = '$Name$'
!-----------------------------------------------------------------------

contains

!#######################################################################

! <SUBROUTINE NAME="atmos_tracer_driver">
!   <OVERVIEW>
!     A routine which allows tracer code to be called.
!   </OVERVIEW>
!   <DESCRIPTION>
!     This subroutine calls the source sink routines for atmospheric
!     tracers. This is the interface between the dynamical core of the
!     model and the tracer code. It should supply all the necessary
!     information to a user that they need in order to calculate the
!     tendency of that tracer with respect to emissions or chemical losses.
!
!   </DESCRIPTION>
!   <TEMPLATE>
!     call atmos_tracer_driver (is, ie, js, je, Time, lon, lat, land, phalf, pfull, r,  &
!                           u, v, t, q, u_star, rdt, rm, rdiag, kbot)
!   </TEMPLATE>
!   <IN NAME="is, ie, js, je" TYPE="integer">
!     Local domain boundaries.
!   </IN>
!   <IN NAME="Time" TYPE="type(time_type)">
!     Model time.
!   </IN>
!   <IN NAME="lon" TYPE="real" DIM="(:,:)">
!     Longitude of the centre of the model gridcells
!   </IN>
!   <IN NAME="lat" TYPE="real" DIM="(:,:)">
!     Latitude of the centre of the model gridcells
!   </IN>
!   <IN NAME="land" TYPE="logical" DIM="(:,:)">
!     fraction of land in grid cell.
!   </IN>
!   <IN NAME="phalf" TYPE="real" DIM="(:,:,:)">
!     Pressures on the model half levels.
!   </IN>
!   <IN NAME="pfull" TYPE="real" DIM="(:,:,:)">
!     Pressures on the model full levels.
!   </IN>
!   <IN NAME="r" TYPE="real" DIM="(:,:,:,:)">
!     The tracer array in the component model.
!   </IN>
!   <IN NAME="u" TYPE="real" DIM="(:,:,:)">
!     Zonal wind speed.
!   </IN>
!   <IN NAME="v" TYPE="real" DIM="(:,:,:)">
!     Meridonal wind speed.
!   </IN>
!   <IN NAME="t" TYPE="real" DIM="(:,:,:)">
!     Temperature.
!   </IN>
!   <IN NAME="q" TYPE="real" DIM="(:,:,:)">
!     Specific humidity. This may also be accessible as a
!                        portion of the tracer array.
!   </IN>
!   <IN NAME="u_star" TYPE="real" DIM="(:,:)">
!     Friction velocity ::
!     The magnitude of the wind stress is density*(ustar**2)
!     The drag coefficient for momentum is u_star**2/(u**2+v**2)
!   </IN>
!   <INOUT NAME="rdt" TYPE="real" DIM="(:,:,:,:)">
!     The tendency of the tracer array in the compenent
!     model. The tendency due to sources and sinks computed
!     in the individual tracer routines should be added to
!     this array before exiting tracer_driver.
!   </INOUT>
!   <IN NAME="rm" TYPE="real" DIM="(:,:,:,:)">
!     The tracer array in the component model for the previous timestep.
!   </IN>
!++amf
!   <IN NAME="flux_sw_down_vis_dir" TYPE="real" DIM="(:,:)">
!     Visible direct radiation at the surface in W / m2
!   </IN>
!   <IN NAME="flux_sw_down_vis_dif" TYPE="real" DIM="(:,:)">
!     Visible diffuse radiation at the surface in W / m2
!   </IN>
!--amf
!   <INOUT NAME="rdiag" TYPE="real" DIM="(:,:,:,:)">
!     The array of diagnostic tracers. As these may be changed within the
!     tracer routines for diagnostic purposes, they need to be writable.
!   </INOUT>
!   <IN NAME="kbot" TYPE="integer, optional" DIM="(:,:)">
!     Integer array describing which model layer intercepts the surface.
!   </IN>
 subroutine atmos_tracer_driver (is, ie, js, je, Time, lon, lat,  &
                           area, z_pbl, rough_mom, &
                           frac_open_sea,land,&
                           phalf, pfull,           &
                           u, v, t, q, r,          &
                           rm, rdt, rdiag, dt,            &
                           u_star, b_star, q_star, &
                           z_half, z_full,         &
                           t_surf_rad, albedo,     &
                           Time_next,              &
                           flux_sw_down_vis_dir,   &
                           flux_sw_down_vis_dif,   &
                           mask,                   &
                           kbot)

!-----------------------------------------------------------------------
integer, intent(in)                           :: is, ie, js, je
type(time_type), intent(in)                   :: Time
real, intent(in),    dimension(:,:)           :: lon, lat
real, intent(in),    dimension(:,:)           :: u_star, b_star, q_star
real, intent(in),    dimension(:,:)           :: land
real, intent(in),    dimension(:,:)           :: area, z_pbl, rough_mom
real, intent(in),    dimension(:,:)           :: frac_open_sea
real, intent(in),    dimension(:,:,:)         :: phalf, pfull
real, intent(in),    dimension(:,:,:)         :: u, v, t, q
real, intent(in),    dimension(:,:,:,:)       :: r
real, intent(in),    dimension(:,:,:,:)       :: rm
real, intent(inout), dimension(:,:,:,:)       :: rdt
real, intent(inout), dimension(:,:,:,ntp+1:)  :: rdiag
real, intent(in)                              :: dt !timestep(used in chem_interface)
real, intent(in),    dimension(:,:,:)         :: z_half !height in meters at half levels
real, intent(in),    dimension(:,:,:)         :: z_full !height in meters at full levels
real, intent(in),    dimension(:,:)           :: t_surf_rad !surface temperature
real, intent(in),    dimension(:,:)           :: albedo
real, intent(in), dimension(:,:)              :: flux_sw_down_vis_dir
real, intent(in), dimension(:,:)              :: flux_sw_down_vis_dif
type(time_type), intent(in)                   :: Time_next
integer, intent(in), dimension(:,:), optional :: kbot
real, intent(in), dimension(:,:,:),  optional :: mask

!-----------------------------------------------------------------------
! Local variables
!-----------------------------------------------------------------------
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtnd, pwt, ozone, o3_prod, &
                                                  aerosol, rho
real, dimension(size(r,1),size(r,2),size(r,3),5) :: rt_sulfate, tr_sulfate
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtndso2, rtndso4,rtnddms
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtndbcphob, rtndbcphil
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtndomphob, rtndomphil
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtndco2, rtndco2_emis
real, dimension(size(r,1),size(r,2),size(rdt,4)) :: dsinku
real, dimension(size(r,1),size(r,2)) ::  w10m_ocean, w10m_land
integer :: year,month,day,hour,minute,second
integer :: jday
real, dimension(size(rdt,1),size(rdt,2),size(rdt,3),size(rdt,4)) :: chem_tend
real, dimension(size(r,1),size(r,2))           :: coszen, fracday, half_day
real :: rrsun
real, dimension(size(r,1),size(r,2),size(r,3)) :: cldf ! cloud fraction
real, dimension(size(r,1),size(r,2),size(r,3)) :: rh  ! relative humidity
real, dimension(size(r,1),size(r,2),size(r,3)) :: lwc ! liq water content
real, dimension(size(r,1),size(r,2),size(r,3)) :: fliq! liq/lwc (f1p)
real, dimension(size(r,1),size(r,2),size(r,3),nt) :: tracer, tracer_orig
real, dimension(size(r,1),size(r,3)) :: dp, temp
real, dimension(size(r,1),size(r,2)) ::  all_salt_settl, all_dust_settl
real, dimension(size(r,1),size(r,2)) ::  suma, ocn_flx_fraction, sum_n_ddep, sum_n_red_ddep, sum_n_ox_ddep
real, dimension(size(r,1),size(r,2)) ::  frland, frsnow, frsea, frice
real, dimension(size(r,1),size(r,2),size(r,3)) :: sumb
integer, dimension(size(r,1),size(r,2)) ::  tropopause_ind

real, dimension(size(r,1),size(r,2),size(r,3)) :: PM1, PM25, PM10

integer :: isulf, i, j, k, id, jd, kd, ntcheck
integer :: nqq  ! index of specific humidity
integer :: nql  ! index of cloud liquid specific humidity
integer :: nqi  ! index of cloud ice water specific humidity
integer :: nqa  ! index of cloud amount
integer :: n, nnn
logical :: used
integer, dimension(6) :: itime   ! JA's time (simpler than model time)

character(len=32) :: tracer_units, tracer_name
real    :: gmt,local_angle
integer :: hh
real :: local_hour_3d(size(r,1),size(r,2),size(r,3)),local_hour
logical :: mask_local_hour(size(r,1),size(r,2),size(r,3))


!-----------------------------------------------------------------------

!   <ERROR MSG="tracer_driver_init must be called first." STATUS="FATAL">
!     Tracer_driver_init needs to be called before tracer_driver.
!   </ERROR>
      if (.not. module_is_initialized)  &
      call error_mesg ('Tracer_driver','tracer_driver_init must be called first.', FATAL)

!-----------------------------------------------------------------------
     id=size(r,1);jd=size(r,2);kd=size(r,3); ntcheck=size(r,4)+size(rdiag,4)
     if (nt /= ntcheck) &
         call error_mesg ('Tracer_driver','number of input tracers not equal total tracers.', FATAL)

      nqa = get_tracer_index(MODEL_ATMOS,'cld_amt')
      nqi = get_tracer_index(MODEL_ATMOS,'ice_wat')
      nql = get_tracer_index(MODEL_ATMOS,'liq_wat')
      nqq = get_tracer_index(MODEL_ATMOS,'sphum')

      rt_sulfate(:,:,:,:)=0.0
      tr_sulfate(:,:,:,:)=0.0
!------------------------------------------------------------------------
! Make local copies of all the tracers
!------------------------------------------------------------------------
!++lwh
      if (use_tau) then
        do n = 1,nt
          if (n <= ntp) then
            tracer(:,:,:,n) = r(:,:,:,n)
          else
            tracer(:,:,:,n) = rdiag(:,:,:,n)
          endif
!------------------------------------------------------------------------
! For tracers other than specific humdity, cloud amount, ice water and &
! liquid water fill eventual negative values
!------------------------------------------------------------------------
!        does tracer need to be adjusted to remain positive definite?
!         if (n /= nqq .and. n/=nqa .and. n/=nqi .and. n/=nql) then
          if ( adjust_positive_def(MODEL_ATMOS,n) ) then
!
            do j=1,jd
              do k=1,kd
                temp(:,k) = tracer(:,j,k,n)
                dp(:,k) = phalf(:,j,k+1)-phalf(:,j,k)
              enddo
              call sjl_fillz(id,kd,1,temp,dp)
              do k=1,kd
                tracer(:,j,k,n) = temp(:,k)
              enddo
            enddo
            if (n <= ntp) then
               rdt(:,:,:,n) = rdt(:,:,:,n) + (tracer(:,:,:,n)-r(:,:,:,n))/dt
            end if
          endif
        end do
      else
        do n = 1,nt
          if (n <= ntp) then
             tracer(:,:,:,n)=rm(:,:,:,n)+rdt(:,:,:,n)*dt
          else
             tracer(:,:,:,n)=rdiag(:,:,:,n)
          end if
!        does tracer need to be adjusted to remain positive definite?
!         if (n /= nqq .and. n/=nqa .and. n/=nqi .and. n/=nql) then
          if ( adjust_positive_def(MODEL_ATMOS,n) ) then
            do j=1,jd
              do k=1,kd
                temp(:,k) = tracer(:,j,k,n)
                dp(:,k) = phalf(:,j,k+1)-phalf(:,j,k)
              enddo
              call sjl_fillz(id,kd,1,temp,dp)
              do k=1,kd
                tracer(:,j,k,n) = temp(:,k)
              enddo
            enddo
            if (n <= ntp) then
               rdt(:,:,:,n) = (tracer(:,:,:,n) - rm(:,:,:,n)) /dt
            end if
          end if
        end do
      end if
!--lwh

      tracer_orig = tracer

!------------------------------------------------------------------------
! Rediagnose meteoroligical variables. Note these parameterizations
! are not consistent with those used elsewhere in the GCM
!------------------------------------------------------------------------

!-----------------------------------------------------------------------
!-----------Calculate relative humidity
!-----------------------------------------------------------------------
      call get_rh(t,q,pfull,rh,mask)
!-----------------------------------------------------------------------
!--------- Calculate wind speed at 10 meters
!-----------------------------------------------------------------------
      call get_w10m(z_full(:,:,kd) - z_half(:,:,kd+1), &
                    u(:,:,kd), v(:,:,kd), &
                    rough_mom, u_star, b_star, q_star, &
                    w10m_ocean, w10m_land, Time_next, is, js)
!-----------------------------------------------------------------------
!------Cloud liquid water content
!-----------------------------------------------------------------------
      if (nqi > 0) then
        lwc(:,:,:)=max(tracer(:,:,:,nqi),0.)
      else
        lwc(:,:,:) = 0.0
      endif
      if (nql > 0) lwc(:,:,:) = lwc(:,:,:) + max(tracer(:,:,:,nql),0.)
      fliq = max(tracer(:,:,:,nql),0.)/max(lwc(:,:,:),1.e-10)
!-----------------------------------------------------------------------
!--------- Cloud fraction -----------------------------
!-----------------------------------------------------------------------
      if (nqa > 0 ) then
!-----  cloud fraction is a prognostic variable ------------
        cldf(:,:,:)= max(0.,min(1.,tracer(:,:,:,nqa) ))
      else
!-----   cloud fraction estimated from RH-------------
        call get_cldf(phalf(:,:,kd+1),pfull,rh,cldf)
      endif
!-----------------------------------------------------------------------
!--------- Get Julian data
!-----------------------------------------------------------------------
      call get_date_julian(Time, year, month, jday, hour, minute, second)
!-----------------------------------------------------------------------
!--------- Get current date
!-----------------------------------------------------------------------
      call get_date(Time, year, month, day, hour, minute, second)

! Calculate cosine of solar zenith angle
!
       call diurnal_solar( lat, lon, Time, coszen, fracday, &
                           rrsun, dt_time=real_to_time_type(dt), &
                           half_day_out=half_day )

!------------------------------------------------------------------------
! Get air mass in layer (in kg/m2), equal to dP/g
!------------------------------------------------------------------------
      do k=1,kd
         pwt(:,:,k)=(phalf(:,:,k+1)-phalf(:,:,k))/grav
         rho(:,:,k) = pwt(:,:,k)/(z_half(:,:,k) - z_half(:,:,k+1))
      enddo

!------------------------------------------------------------------------
!    define various land fractions needed for dry deposition calculation.
!------------------------------------------------------------------------
      frland(:,:) = min(1., max(0.,     land(:,:) ) )
      frsea(:,:)  = min(1., max(0., 1. - frland(:,:) ) )
!      frsnow(:,:) = min(frland(:,:), max(0., snow_area(:,:) ) )
      frice(:,:)  = min(frsea(:,:),   &
                               max(0., frsea(:,:) - frac_open_sea(:,:) ) )

!------------------------------------------------------------------------
!    output land fraction information, if desired.
!------------------------------------------------------------------------
   used = send_data ( id_landfr, frland, Time_next, is_in =is,js_in=js)
   used = send_data ( id_seaicefr, frice, Time_next, is_in =is,js_in=js)
!   used = send_data ( id_snowfr, frsnow, Time_next, is_in =is,js_in=js)
!   used = send_data ( id_vegnfr, vegn_cover, Time_next, is_in =is,js_in=js)
!   used = send_data ( id_vegnlai, vegn_lai, Time_next, is_in =is,js_in=js)

!------------------------------------------------------------------------
! For tracers other than specific humdity, cloud amount, ice water and &
! liquid water calculate flux at surface due to dry deposition
!------------------------------------------------------------------------
!++lwh

     sum_n_ddep(:,:) = 0.
     sum_n_ox_ddep(:,:) = 0.
     sum_n_red_ddep(:,:) = 0.

      do n=1,ntp
         if (n /= nqq .and. n/=nqa .and. n/=nqi .and. n/=nql) then
            call dry_deposition( n, is, js, u(:,:,kd), v(:,:,kd), t(:,:,kd), &
                                 pwt(:,:,kd), pfull(:,:,kd), &
                                 z_half(:,:,kd)-z_half(:,:,kd+1), u_star, &
                                 land, dsinku(:,:,n), dt, &
                                 tracer(:,:,kd,n), Time, Time_next, &
                                 lon, half_day, &
                                 drydep_data(n))!, frland, frice, frsnow, &
!                                 vegn_cover, vegn_lai, &
!                                 b_star, z_pbl, rough_mom)
            rdt(:,:,kd,n) = rdt(:,:,kd,n) - dsinku(:,:,n)
            if ( step_update_tracer ) then
               tracer(:,:,kd,n) = tracer(:,:,kd,n) - dsinku(:,:,n)*dt
            end if
            if (id_n_ddep>0) sum_n_ddep = sum_n_ddep + pwt(:,:,kd)*dsinku(:,:,n)*14./wtmair*nb_n(n)
            if (id_n_ox_ddep>0) sum_n_ox_ddep = sum_n_ox_ddep + pwt(:,:,kd)*dsinku(:,:,n)*14./wtmair*nb_n_ox(n)
            if (id_n_red_ddep>0) sum_n_red_ddep = sum_n_red_ddep + pwt(:,:,kd)*dsinku(:,:,n)*14./wtmair*nb_n_red(n)

            if (id_tracer_ddep_kg_m2_s(n)>0) then
               used = send_data ( id_tracer_ddep_kg_m2_s(n), dsinku(:,:,n)*pwt(:,:,kd)*conv_vmr_mmr(n), Time_next, is_in=is,js_in=js)
            end if

         end if
      enddo

      if (id_om_ddep > 0 .and. nomphilic > 0 .and. nomphobic > 0) then
        used  = send_data (id_om_ddep,  &
         pwt(:,:,kd)*(dsinku(:,:,nomphilic) + dsinku(:,:,nomphobic)),  &
                                              Time_next, is_in=is, js_in=js)
      endif
      if (id_bc_ddep > 0 .and. nbcphilic > 0 .and. nbcphobic > 0) then
        used  = send_data (id_bc_ddep,  &
         pwt(:,:,kd)*(dsinku(:,:,nbcphilic) + dsinku(:,:,nbcphobic)),  &
                                               Time_next, is_in=is, js_in=js)
      endif
      if (id_nh4_ddep_cmip > 0 .and. nNH4NO3 > 0 .and. nNH4 > 0) then
        used  = send_data (id_nh4_ddep_cmip,  &
        0.018*1.0e03*pwt(:,:,kd)*(dsinku(:,:,nNH4NO3) + dsinku(:,:,nNH4))/WTMAIR,  &
                                              Time_next, is_in=is, js_in=js)
      endif


      !---- cmip variables ----
      if (id_n_ox_ddep > 0) used = send_data (id_n_ox_ddep, sum_n_ox_ddep, Time_next, &
                                              is_in=is, js_in=js)

      if (id_dryso2 > 0) then
        if (nSO2_cmip > 0) then
          used = send_data (id_dryso2, 0.064*1.e3*pwt(:,:,kd)*dsinku(:,:,nSO2_cmip)/WTMAIR, &
                          Time_next, is_in=is, js_in=js)
        else if (nSO2 > 0) then ! fast-aerosol simpleSO2
          used = send_data (id_dryso2, 0.064*1.e3*pwt(:,:,kd)*dsinku(:,:,nSO2)/WTMAIR, &
                          Time_next, is_in=is, js_in=js)
        endif
      endif
      if (id_dryso4 > 0) then
        if (nSO4_cmip > 0) then
          used = send_data (id_dryso4, 0.096*1.e3*pwt(:,:,kd)*dsinku(:,:,nSO4_cmip)/WTMAIR, &
                            Time_next, is_in=is, js_in=js)
        else if (nSO4 > 0) then ! fast-aerosol simpleSO4
          used = send_data (id_dryso4, 0.096*1.e3*pwt(:,:,kd)*dsinku(:,:,nSO4)/WTMAIR, &
                            Time_next, is_in=is, js_in=js)
        endif
      endif
      if (id_drydms > 0) then
        if (nDMS_cmip > 0) then
          used = send_data (id_drydms, 0.062*1.e3*pwt(:,:,kd)*dsinku(:,:,nDMS_cmip)/WTMAIR, &
                          Time_next, is_in=is, js_in=js)
        else if (nDMS > 0) then ! fast-aerosol simpleDMS
          used = send_data (id_drydms, 0.062*1.e3*pwt(:,:,kd)*dsinku(:,:,nDMS)/WTMAIR, &
                          Time_next, is_in=is, js_in=js)
        endif
      endif
      if (id_drynh3 > 0 .and. nNH3_cmip > 0) then
        used = send_data (id_drynh3, 0.017*1.e3*pwt(:,:,kd)*dsinku(:,:,nNH3_cmip)/WTMAIR, &
                          Time_next, is_in=is, js_in=js)
      endif
      if (id_drysoa > 0 .and. nSOA > 0) then
        used = send_data (id_drysoa, pwt(:,:,kd)*dsinku(:,:,nSOA), Time_next, is_in=is, js_in=js)
      endif

      if (id_drynh4 > 0 .and. nNH4NO3 > 0 .and. nNH4 > 0) then  ! same as nh4_ddep_cmip
        used  = send_data (id_drynh4,  &
                   0.018*1.e3*pwt(:,:,kd)*(dsinku(:,:,nNH4NO3)+dsinku(:,:,nNH4))/WTMAIR,  &
                              Time_next, is_in=is, js_in=js)
      endif
      if (id_drybc > 0 .and. nbcphilic > 0 .and. nbcphobic > 0) then  ! same as bc_ddep
        used  = send_data (id_drybc,  &
           pwt(:,:,kd)*(dsinku(:,:,nbcphilic) + dsinku(:,:,nbcphobic)),  &
                                               Time_next, is_in=is, js_in=js)
      endif
      if (id_drypoa > 0 .and. nomphilic > 0 .and. nomphobic > 0) then
        used  = send_data (id_drypoa,  &
            pwt(:,:,kd)*(dsinku(:,:,nomphilic) + dsinku(:,:,nomphobic)),  &
                                     Time_next, is_in=is, js_in=js)
      endif
      if (id_dryoa > 0 .and. nomphilic > 0 .and. nomphobic > 0 .and. nSOA > 0) then
        used  = send_data (id_dryoa,  &
            pwt(:,:,kd)*(dsinku(:,:,nomphilic) + dsinku(:,:,nomphobic) + dsinku(:,:,nSOA)),  &
                                     Time_next, is_in=is, js_in=js)
      endif

      do n=1,ntp
         if (id_tracer_col_kg_m2(n).gt.0) then
            suma = 0.
            do k=1,kd
               suma(:,:) = suma(:,:) + pwt(:,:,k)*tracer(:,:,k,n)
            end do
            suma(:,:) = conv_vmr_mmr(n)*suma(:,:)
            used      = send_data (id_tracer_col_kg_m2(n), suma, Time_next, is_in=is, js_in=js)
         end if
      end do

      if (id_bc_col_kg_m2.gt.0 .and. nomphilic.gt.0 .and. nomphobic.gt.0) then
        suma = 0.
        do k=1,kd
           suma(:,:) = suma(:,:) + pwt(:,:,k)*(tracer(:,:,k,nbcphilic)+tracer(:,:,k,nbcphobic))
        end do
        used = send_data (id_bc_col_kg_m2, suma, Time_next, is_in=is, js_in=js)
      end if

      if (id_dust_col_kg_m2.gt.0) then
        suma = 0.
        do n=1,ntp
          if (is_dust_tracer(n)) then
             do k=1,kd
                suma(:,:) = suma(:,:) + pwt(:,:,k)*tracer(:,:,k,n)
             end do
          end if
        end do
        used = send_data (id_dust_col_kg_m2, suma, Time_next, is_in=is, js_in=js)
      end if

      if (id_seasalt_col_kg_m2.gt.0) then
        suma = 0.
        do n=1,ntp
          if (is_seasalt_tracer(n)) then
             do k=1,kd
                suma(:,:) = suma(:,:) + pwt(:,:,k)*tracer(:,:,k,n)
             end do
          end if
        end do
        used = send_data (id_seasalt_col_kg_m2, suma, Time_next, is_in=is, js_in=js)
      end if

!----------------------------------------------------------------------
!   output the nh4no3 and nh4 loads.
!----------------------------------------------------------------------
!THIS IS WRRONG IN AM4
      if (nNH4NO3 > 0 .and. nNH4 > 0) then
        if (id_nh4_col > 0 .or. id_loadnh4 > 0) then
          suma = 0.
          do k=1,kd
            suma(:,:) = suma(:,:) + pwt(:,:,k)*(tracer(:,:,k,nNH4) + &
                             tracer(:,:,k,nNH4NO3))
          end do
          suma(:,:) = 0.018*1.0e03*suma(:,:)/WTMAIR
          if (id_nh4_col > 0) then
            used  = send_data (id_nh4_col, suma, Time_next, is_in=is, js_in=js)
          endif
          ! cmip named field
          if (id_loadnh4 > 0) then
            used  = send_data (id_loadnh4, suma, Time_next, is_in=is, js_in=js)
          endif
        endif
      endif

      if (nNH4NO3 > 0) then
        if (id_nh4no3_col > 0 .or. id_loadno3 > 0) then
          suma = 0.
          do k=1,kd
            suma(:,:) = suma(:,:) + pwt(:,:,k)*tracer(:,:,k,nNH4NO3)
          end do
          suma(:,:) = 0.062*1.0e03*suma(:,:)/WTMAIR
          if (id_nh4no3_col > 0) then
            used  = send_data (id_nh4no3_col, suma(:,:), Time_next, is_in=is, js_in=js)
          endif
          ! cmip named field
          if (id_loadno3 > 0) then
            used  = send_data (id_loadno3, suma(:,:), Time_next, is_in=is, js_in=js)
          endif
        endif
      endif

!----------------------------------------------------------------------
!   output the tracer fields needed for CMIP.
!----------------------------------------------------------------------
      if (id_nh4_cmip > 0) then
        used  = send_data (id_nh4_cmip,  &
               0.018*1.0e03* (tracer(:,:,:,nNH4NO3) + &
                              tracer(:,:,:,nNH4)) /WTMAIR,  &
                                          Time_next, is_in=is, js_in=js, ks_in=1)
      endif
      if (id_nh4_cmipv2 > 0) then
        used  = send_data (id_nh4_cmipv2,  &
               0.018*1.0e03*rho(:,:,:)* (tracer(:,:,:,nNH4NO3) + &
                              tracer(:,:,:,nNH4)) /WTMAIR,  &
                                          Time_next, is_in=is, js_in=js, ks_in=1)
      endif
      if(id_nh4no3_cmip > 0) then
        used  = send_data (id_nh4no3_cmip,  &
                0.062*1.0e03*tracer(:,:,:,nNH4NO3)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
      if(id_nh4no3_cmipv2 > 0) then
        used  = send_data (id_nh4no3_cmipv2,  &
                0.062*1.0e03*rho(:,:,:)*tracer(:,:,:,nNH4NO3)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
     if(id_so2_cmip > 0) then
       used  = send_data (id_so2_cmip,  &
                 0.064*1.0e03*tracer(:,:,:,nSO2_cmip)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
     if(id_so2_cmipv2 > 0) then
       used  = send_data (id_so2_cmipv2,  &
               0.064*1.0e03*rho(:,:,:)*tracer(:,:,:,nSO2_cmip)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
     if(id_dms_cmip > 0) then
       used  = send_data (id_dms_cmip,  &
                0.062*1.0e03*tracer(:,:,:,nDMS_cmip)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
     if(id_dms_cmipv2 > 0) then
       used  = send_data (id_dms_cmipv2,  &
                0.062*1.0e03*rho(:,:,:)*tracer(:,:,:,nDMS_cmip)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif

     !---- cmip named variables ----
     if (query_cmip_diag_id(ID_airmass)) then
         used = send_cmip_data_3d (ID_airmass,  &
                 pwt, Time_next, is_in=is, js_in=js, ks_in=1)
     endif




     if (nNH4NO3 > 0 .and. nNH4 > 0) then
       ! concentration
       if (query_cmip_diag_id(ID_concnh4)) then
         used = send_cmip_data_3d (ID_concnh4,  &
                 0.018*1.0e03*rho(:,:,:)*(tracer(:,:,:,nNH4NO3) + tracer(:,:,:,nNH4)) /WTMAIR,  &
                                   Time_next, is_in=is, js_in=js, ks_in=1)
       endif
       ! surface concentration (lowest level)
       if (id_sconcnh4 > 0) then
         used = send_data (id_sconcnh4, &
             0.018*1.0e03*rho(:,:,kd)*(tracer(:,:,kd,nNH4NO3) + tracer(:,:,kd,nNH4)) /WTMAIR,  &
                                 Time_next, is_in=is, js_in=js)
       endif
     endif

     if (nNH4NO3 > 0) then
       if (query_cmip_diag_id(ID_concno3)) then
         used = send_cmip_data_3d ( ID_concno3,  &
                0.062*1.0e03*rho(:,:,:)*tracer(:,:,:,nNH4NO3)/WTMAIR,  &
                                 Time_next, is_in=is, js_in=js, ks_in=1)
       endif
       ! cmip surface concentration (lowest level)
       if (id_sconcno3 > 0) then
         used = send_data (id_sconcno3, &
                0.062*1.0e03*rho(:,:,kd)*tracer(:,:,kd,nNH4NO3)/WTMAIR,  &
                                 Time_next, is_in=is, js_in=js)
       endif
     endif

     if (query_cmip_diag_id(ID_concso2)) then
       if (nSO2_cmip > 0) then
         used = send_cmip_data_3d ( ID_concso2, tracer(:,:,:,nSO2_cmip), &
                                    Time_next, is_in=is, js_in=js, ks_in=1)
       else if (nSO2 > 0) then ! fast-aerosol simpleSO2
         used = send_cmip_data_3d ( ID_concso2, tracer(:,:,:,nSO2), &
                                    Time_next, is_in=is, js_in=js, ks_in=1)
       endif
     endif
     if (query_cmip_diag_id(ID_concdms)) then
       if (nDMS_cmip > 0) then
         used = send_cmip_data_3d ( ID_concdms, tracer(:,:,:,nDMS_cmip), &
                                Time_next, is_in=is, js_in=js, ks_in=1)
       else if (nDMS > 0) then ! fast-aerosol simpleDMS
         used = send_cmip_data_3d ( ID_concdms, tracer(:,:,:,nDMS), &
                                Time_next, is_in=is, js_in=js, ks_in=1)
       endif
     endif

     PM25 = 0.
     PM1  = 0.
     PM10 = 0.


     do n=1,nt

        if (frac_pm25(n).gt.0.)  PM25 = PM25+conv_vmr_mmr(n)*tracer(:,:,:,n)*frac_pm25(n)
        if (frac_pm1(n).gt.0.)   PM1  = PM1+conv_vmr_mmr(n)*tracer(:,:,:,n)*frac_pm1(n)
        if (frac_pm10(n).gt.0.)  PM10 = PM10+conv_vmr_mmr(n)*tracer(:,:,:,n)*frac_pm10(n)

        if ( query_cmip_diag_id(ID_tracer_mol_mol(n)) ) then
           used = send_cmip_data_3d ( ID_tracer_mol_mol(n), tracer(:,:,:,n), &
                Time_next, is_in=is, js_in=js, ks_in=1)
        end if
        if ( id_tracer_surf_mol_mol(n) .gt. 0 ) then
           used = send_data ( id_tracer_surf_mol_mol(n), tracer(:,:,kd,n), &
                Time_next, is_in=is, js_in=js)
        end if
        if ( query_cmip_diag_id(ID_tracer_kg_kg(n)) ) then
           used = send_cmip_data_3d ( ID_tracer_kg_kg(n), conv_vmr_mmr(n)*tracer(:,:,:,n), &
                Time_next, is_in=is, js_in=js, ks_in=1)
        end if
        if ( id_tracer_surf_kg_kg(n) .gt. 0 ) then
           used = send_data ( id_tracer_surf_kg_kg(n), conv_vmr_mmr(n)*tracer(:,:,kd,n), &
                Time_next, is_in=is, js_in=js)
        end if

     end do

     if (id_pm25_surf .gt. 0) then
        used = send_data ( id_pm25_surf, PM25(:,:,kd), &
             Time_next, is_in=is, js_in=js)
     end if

     if ( query_cmip_diag_id(ID_OM) .and. nomphilic > 0 .and. nomphobic > 0) then
        used = send_cmip_data_3d ( ID_OM, tracer(:,:,:,nomphilic)+tracer(:,:,:,nomphobic), &
             Time_next, is_in=is, js_in=js, ks_in=1)
     end if

     if ( query_cmip_diag_id(ID_BC) .and. nbcphilic > 0 .and. nbcphobic > 0) then
        used = send_cmip_data_3d ( ID_BC, tracer(:,:,:,nbcphilic)+tracer(:,:,:,nbcphobic), &
             Time_next, is_in=is, js_in=js, ks_in=1)
     end if

     if ( query_cmip_diag_id(ID_DUST) ) then
        sumb = 0.
        do n=1,ntp
          if (is_dust_tracer(n)) then
             sumb(:,:,:) = sumb(:,:,:) + tracer(:,:,:,n)
          end if
        end do
        used = send_cmip_data_3d ( ID_DUST, sumb(:,:,:), &
             Time_next, is_in=is, js_in=js, ks_in=1)
     end if

     if ( query_cmip_diag_id(ID_SS) ) then
        sumb = 0.
        do n=1,ntp
          if (is_seasalt_tracer(n)) then
             sumb(:,:,:) = sumb(:,:,:) + tracer(:,:,:,n)
          end if
        end do
        used = send_cmip_data_3d ( ID_SS, sumb(:,:,:), &
             Time_next, is_in=is, js_in=js, ks_in=1)
     end if

     if ( query_cmip_diag_id(ID_pm25)) then
        used = send_cmip_data_3d ( ID_pm25, pm25, &
             Time_next, is_in=is, js_in=js, ks_in=1)
     end if

     if ( query_cmip_diag_id(ID_pm1)) then
        used = send_cmip_data_3d ( ID_pm1, pm1, &
             Time_next, is_in=is, js_in=js, ks_in=1)
     end if

     if ( query_cmip_diag_id(ID_pm10)) then
        used = send_cmip_data_3d ( ID_pm10, pm10, &
             Time_next, is_in=is, js_in=js, ks_in=1)
     end if


!------------------------------------------------------------------------
! Compute tropopause diagnostics
!------------------------------------------------------------------------
   call mpp_clock_begin (tropopause_clock)
   call atmos_tropopause(is, ie, js, je, Time, Time_next, t, pfull, z_full,     &
                         tropopause_ind)
   call mpp_clock_end (tropopause_clock)

      if (id_toz > 0 .and. no3 > 0) then
        suma = 0.
        do k=1,kd
           suma(:,:) = suma(:,:) + pwt(:,:,k)*tracer(:,:,k,no3)
        end do
	suma = suma * AVOGNO / (WTMAIR*1.e-3 * o3_column_factor)
        used = send_data (id_toz, suma, Time_next, is_in=is, js_in=js)
      end if

      if (id_tropoz > 0 .and. no3 > 0) then
        suma = 0.
	do i = 1,id
	do j = 1,jd
           do k=1,tropopause_ind(i,j)
              suma(:,:) = suma(:,:) + pwt(:,:,k)*tracer(:,:,k,no3)
           end do
        end do
        end do
	suma = suma * AVOGNO / (WTMAIR*1.e-3 * o3_column_factor)
        used = send_data (id_tropoz, suma, Time_next, is_in=is, js_in=js)
      end if

!------------------------------------------------------------------------
! Compute radon source-sink tendency
!------------------------------------------------------------------------
    call mpp_clock_begin (radon_clock)
    do nnn = 1, size(nradon(:))
     if (nradon(nnn) > 0) then
       if (nradon(nnn) > nt) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for radon', FATAL)
         call atmos_radon_sourcesink (lon,lat,land,pwt,tracer(:,:,:,nradon(nnn)),  &
                                 rtnd, Time, kbot)
       rdt(:,:,:,nradon(nnn))=rdt(:,:,:,nradon(nnn))+rtnd(:,:,:)
       if ( step_update_tracer ) then
          tracer(:,:,:,nradon(nnn)) = tracer(:,:,:,nradon(nnn)) + rtnd(:,:,:)*dt
       end if
    endif

   end do
   call mpp_clock_end (radon_clock)

!------------------------------------------------------------------------
! Compute convection tracer source-sink tendency
!------------------------------------------------------------------------
   call mpp_clock_begin (convect_clock)
   do nnn = 1, size(nconvect(:))
     if (nconvect(nnn) > 0) then
       if (nconvect(nnn) > nt) call error_mesg ('Tracer_driver', &
           'Number of tracers .lt. number for convection tracer', FATAL)
       call atmos_cnvct_tracer_sourcesink (lon,lat,land,pwt,  &
                                       tracer(:,:,:,nconvect(nnn)),  &
                                       rtnd, Time, is, ie, js, je,kbot)
       rdt(:,:,:,nconvect(nnn))=rdt(:,:,:,nconvect(nnn))+rtnd(:,:,:)
     endif
   end do
   call mpp_clock_end (convect_clock)

!------------------------------------------------------------------------
! Stratospheric chemistry
!------------------------------------------------------------------------
  if(do_coupled_stratozone) then
    call mpp_clock_begin (stratozone_clock)
    itime(:) = 0
    call get_date(time,itime(1),itime(2),itime(3),itime(4),itime(5),itime(6))

    call strat_chem(lon,lat,r,rdt,pfull,t,itime,is,ie,js,je,dt,coszen,  &
!    ozon,cosp,cosphc,photo,solardata,chlb,ozb,dfdage,tropc,anoy,        &
    nsphum,chem_tend,ozone,o3_prod,aerosol,mpp_pe())
!
! The water vapour tendency is included here, partially coupling the
! chemistry to the radiation
!
      rdt(:,:,:,:) = rdt(:,:,:,:) + chem_tend(:,:,:,1:ntp)
      if (step_update_tracer) then
         tracer(:,:,:,1:ntp) = tracer(:,:,:,1:ntp) + chem_tend(:,:,:,1:ntp)*dt
      end if

      if(nt.gt.(ntp+1))  then
! Modify the diagnostic tracers.
        tracer(:,:,:,no3)      = ozone(:,:,:)
        tracer(:,:,:,no3ch)    = o3_prod(:,:,:)
        if (naerosol > 0 ) tracer(:,:,:,naerosol) = aerosol(:,:,:)
      endif
    call mpp_clock_end (stratozone_clock)
  endif

!-----------------------------------------------------------------------
!   define sfc area through which tracer flux from ocean is allowed; it
!   is area fraction of sea unless prevent_flux_through_ice is .true.,
!   in which case it is the fractional area of open sea.
!-----------------------------------------------------------------------
   if (prevent_flux_through_ice) then
     ocn_flx_fraction = frac_open_sea
   else
     ocn_flx_fraction = 1. - land
   endif

!------------------------------------------------------------------------
! Tropospheric chemistry
!------------------------------------------------------------------------
   if ( do_tropchem ) then
!------------------------------------------------------------------------
! Compute age tracer source-sink tendency
!------------------------------------------------------------------------
      if (nage > 0) then
        call mpp_clock_begin (age_tracer_clock)
        if (nage > nt) call error_mesg ('Tracer_driver', &
           'Number of tracers .lt. number for age tracer', FATAL)
        call atmos_age_tracer( lon, lat, pwt,  &
                               tracer(:,:,:,nage),  &
                               rtnd, Time, kbot)
          rdt(:,:,:,nage)=rdt(:,:,:,nage)+rtnd(:,:,:)
        call mpp_clock_end (age_tracer_clock)
      endif

      call mpp_clock_begin (tropchem_clock)
      call tropchem_driver( lon, lat, land, ocn_flx_fraction, pwt, &
                            tracer(:,:,:,1:ntp),chem_tend, &
                            Time, phalf, pfull, t, is, ie, js, je, dt, &
                            z_half, z_full, q, t_surf_rad, albedo, coszen, rrsun, &
                            area, w10m_ocean, &
                            flux_sw_down_vis_dir, flux_sw_down_vis_dif, &
                            half_day, &
                            Time_next, tracer(:,:,:,MIN(ntp+1,nt):nt), kbot)
      rdt(:,:,:,:) = rdt(:,:,:,:) + chem_tend(:,:,:,:)
      call mpp_clock_end (tropchem_clock)
   endif

!! RSH 4/8/04
!! note that if there are no diagnostic tracers, that argument in the
!! call to sourcesink should be made optional and omitted in the calls
!! below. note the switch in argument order to make this argument
!! optional.
    if (nt == ntp) then  ! implies no diagnostic tracers
!   if(do_tropchem) then
!      if(present(kbot)) then
!        call sourcesink(lon,lat,land,pwt,r+rdt*dt,chem_tend,Time,phalf,pfull,t,is,js,je,dt,&
!                          z_half, z_full,q,t_surf_rad,albedo,coszen, Time_next,&
!                          u,v,u_star,&
!                          kbot)
!      else
!       call sourcesink(lon,lat,land,pwt,r+rdt*dt,chem_tend,Time,phalf,pfull,t,is,js,je,dt, &
!                          z_half, z_full,q, t_surf_rad, albedo, coszen, Time_next, &
!                          u,v,u_star,&
!                          )
!      endif
!      rdt(:,:,:,:) = rdt(:,:,:,:) + chem_tend(:,:,:,:)
!   endif
      else   ! case of diagnostic tracers being present
!   if(do_tropchem) then
!      if(present(kbot)) then
!        call sourcesink(lon,lat,land,pwt,r+rdt*dt,chem_tend,Time,phalf,pfull,t,is,js,je,dt,&
!                          z_half, z_full,q,t_surf_rad,albedo,coszen, Time_next,&
!                          u,v,u_star,&
!       rdiag=rm(:,:,:,nt+1:ntp),  &  ! (the diagnostic tracers)
!                   kbot=kbot)
!      else
!       call sourcesink(lon,lat,land,pwt,r+rdt*dt,chem_tend,Time,phalf,pfull,t,is,js,je,dt, &
!                          z_half, z_full,q, t_surf_rad, albedo, coszen, Time_next, &
!                          u,v,u_star,&
!        rdiag=rm(:,:,:,nt+1:ntp),  & ! (the diagnostic tracers)
!                          )
!      endif
!      rdt(:,:,:,:) = rdt(:,:,:,:) + chem_tend(:,:,:,:)
!   endif
   endif  ! (no diagnostic tracers)

!------------------------------------------------------------------------
!   carbonaceous aerosols
!------------------------------------------------------------------------
   if (nbcphobic > 0 .and. nbcphilic > 0 .and. &
       nomphobic > 0 .and. nomphilic > 0) then
         if (nbcphobic > ntp .or. nbcphilic > ntp) &
            call error_mesg ('Tracer_driver', &
            'Number of tracers .lt. number for black carbon', FATAL)
     call mpp_clock_begin (carbon_clock)
     call atmos_carbon_aerosol_driver(lon,lat,ocn_flx_fraction, pfull,phalf,z_half,z_pbl, &
                                      t_surf_rad, w10m_ocean, &
                                      T, pwt, &
                                      tracer(:,:,:,nbcphobic), rtndbcphob, &
                                      tracer(:,:,:,nbcphilic), rtndbcphil, &
                                      tracer(:,:,:,nomphobic), rtndomphob, &
                                      tracer(:,:,:,nomphilic), rtndomphil, &
                                      tracer(:,:,:,nOH),    &
                                      Time_next,is,ie,js,je)
      rdt(:,:,:,nbcphobic)=rdt(:,:,:,nbcphobic)+rtndbcphob(:,:,:)
      rdt(:,:,:,nbcphilic)=rdt(:,:,:,nbcphilic)+rtndbcphil(:,:,:)
      rdt(:,:,:,nomphobic)=rdt(:,:,:,nomphobic)+rtndomphob(:,:,:)
      rdt(:,:,:,nomphilic)=rdt(:,:,:,nomphilic)+rtndomphil(:,:,:)
      call mpp_clock_end (carbon_clock)
   endif
!------------------------------------------------------------------------
! Mineral Dust
!------------------------------------------------------------------------
  call mpp_clock_begin (dust_clock)
   if (do_dust) then
      call atmos_dust_sourcesink(lon,lat,land,pwt, dt, &
              z_half, pfull, w10m_land, t, rh, &
              tracer(:,:,:,:), dsinku(:,:,:), rdt(:,:,:,:), &
              Time, is,ie,js,je, kbot)
   endif
   call mpp_clock_end (dust_clock)

!------------------------------------------------------------------------
!sea salt
!------------------------------------------------------------------------
   call mpp_clock_begin (seasalt_clock)
   if (do_seasalt) then
      call atmos_sea_salt_sourcesink(lon,lat,ocn_flx_fraction,pwt, &
              z_half, pfull, w10m_ocean, t, rh, &
              tracer(:,:,:,:), dsinku(:,:,:), rdt(:,:,:,:), dt, &
              Time, is,ie,js,je, kbot)
   endif
   call mpp_clock_end (seasalt_clock)

!------------------------------------------------------------------------
! Sulfur chemistry
!------------------------------------------------------------------------
   if (nDMS > 0 .and. nSO2 > 0 .and. nSO4 > 0) then
      if (nDMS > ntp ) call error_mesg ('Tracer_driver', &
                     'Number of tracers .lt. number for DMS', FATAL)
      if (nSO2 > ntp ) call error_mesg ('Tracer_driver', &
                     'Number of tracers .lt. number for SO2', FATAL)
      if (nSO4 > ntp ) call error_mesg ('Tracer_driver', &
                     'Number of tracers .lt. number for SO4', FATAL)

      call mpp_clock_begin (sulfur_clock)
      if (nDMS > 0 ) then
        call atmos_DMS_emission(lon, lat, area, ocn_flx_fraction, t_surf_rad, &
             w10m_ocean, pwt, rtnddms, Time, Time_next, is,ie,js,je,kbot)
        rdt(:,:,kd,nDMS) = rdt(:,:,kd,nDMS) + rtnddms(:,:,kd)
      endif
      call atmos_SOx_emission(lon, lat, area, land, &
               z_pbl, z_half, phalf, pwt, rtndso2, rtndso4, &
               Time, Time_next, is,ie,js,je,kbot)
      rdt(:,:,:,nSO4) = rdt(:,:,:,nSO4) + rtndso4(:,:,:)
      rdt(:,:,:,nSO2) = rdt(:,:,:,nSO2) + rtndso2(:,:,:)
      do isulf=1,nsulfate
        if (do_tracer_sulfate(isulf)) tr_sulfate(:,:,:,isulf)= &
              tracer(:,:,:,tr_nbr_sulfate(isulf))
      enddo
      call atmos_SOx_chem( pwt, t, pfull, phalf, dt, lwc, fliq, cldf, &
                jday,hour,minute,second,lat,lon,    &
                do_tracer_sulfate, tr_sulfate, rt_sulfate, &
                tracer(:,:,:,nOH), &
                Time,Time_next, is,ie,js,je,kbot)
      do isulf=1,nsulfate
        if (do_tracer_sulfate(isulf)) rdt(:,:,:,tr_nbr_sulfate(isulf))= &
              rdt(:,:,:,tr_nbr_sulfate(isulf))+rt_sulfate(:,:,:,isulf)
      enddo
      call mpp_clock_end (sulfur_clock)
   endif

!------------------------------------------------------------------------
! Secondary organic aerosols
!------------------------------------------------------------------------
   if (nSOA > 0 ) then
      if (nSOA > ntp ) call error_mesg ('Tracer_driver', &
                     'Number of tracers .lt. number for SOA', FATAL)

      call mpp_clock_begin (SOA_clock)
      call atmos_SOA_chem(pwt ,t, pfull, phalf, dt, &
                jday, hour, minute, second, lat, lon,    &
                tracer(:,:,:,nSOA), &
		tracer(:,:,:,nOH), &
		tracer(:,:,:,nC4H10), &
		rtnd, Time, Time_next, is,ie,js,je,kbot )

      rdt(:,:,:,nSOA)=rdt(:,:,:,nSOA)+rtnd(:,:,:)
      call mpp_clock_end (SOA_clock)

   endif

!------------------------------------------------------------------------
! Sulfur hexafluoride (SF6)
!------------------------------------------------------------------------
   if (nsf6 > 0) then
         if (nsf6 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for sulfur hexafluoride', FATAL)
         call mpp_clock_begin (sf6_clock)
         call atmos_sf6_sourcesink (lon,lat,land,pwt,r(:,:,:,nsf6),  &
                                 rtnd, Time,is,ie,js,je, kbot)
      rdt(:,:,:,nsf6)=rdt(:,:,:,nsf6)+rtnd(:,:,:)
         call mpp_clock_end (sf6_clock)
   endif

   if (nch3i > 0) then
         if (nch3i > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers < number for ch3i', FATAL)
         call mpp_clock_begin (ch3i_clock)
         call atmos_ch3i( lon, lat, land, pwt, &
                          r(:,:,:,nch3i)+dt*rdt(:,:,:,nch3i), rtnd,  &
                          Time, phalf, pfull, t, is, js, dt, &
                          z_half, z_full,q, t_surf_rad, albedo, coszen, &
                          Time_next, kbot)
         rdt(:,:,:,nch3i) = rdt(:,:,:,nch3i) + rtnd(:,:,:)
         call mpp_clock_end (ch3i_clock)
   endif

!------------------------------------------------------------------------
! Compute CO2 source-sink tendency:  tracer(:,:,:,nco2) is in moist mmr
!  atmos_co2_sourcesink will convert to dry vmr and return rtnd as moist mmr
!------------------------------------------------------------------------
   if (nco2 > 0 ) then
         if (nco2 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers < number for co2', FATAL)
         call mpp_clock_begin (co2_clock)
         call atmos_co2_emissions (is, ie, js, je, Time, Time_next, dt, pwt, tracer(:,:,:,nco2),     &
                                   tracer(:,:,:,nsphum), rtndco2_emis, kbot)
         rdt(:,:,:,nco2)=rdt(:,:,:,nco2)+rtndco2_emis(:,:,:)

         call atmos_co2_sourcesink (is, ie, js, je, Time, Time_next, dt, pwt, tracer(:,:,:,nco2),     &
                                    tracer(:,:,:,nsphum), rtndco2)
         rdt(:,:,:,nco2)=rdt(:,:,:,nco2)+rtndco2(:,:,:)
         call mpp_clock_end (co2_clock)
   endif

   if (ncodirect > 0 .or. ne90 > 0) then
      if (ncodirect > ntp .or. ne90 > ntp) call error_mesg ('Tracer_driver', &
                         'Number of tracers < number for codirect', FATAL)
      call regional_tracer_driver( lon, lat, pwt, tracer, chem_tend, &
                                   Time, phalf, is, js, kbot)
      rdt(:,:,:,:) = rdt(:,:,:,:) + chem_tend(:,:,:,:)
   endif

!------------------------------------------------------------------------
! Save diagnostic tracer concentrations back to tracer array
!------------------------------------------------------------------------
   if (nt > ntp) then
      rdiag(:,:,:,ntp+1:nt) = tracer(:,:,:,ntp+1:nt)
   end if


   !save tracer diagnostics

!calculate local time
   gmt = universal_time(Time) !time of day midnight = 0
   do j=1,jd
      do i=1,id
         local_angle = gmt + lon(i,j)
         if (local_angle >= 2.*PI) local_angle = local_angle - 2*pi
         local_hour        = local_angle *12./PI
         local_hour_3d(i,j,:) = local_hour
      end do
   end do


   do hh=1,24
      mask_local_hour = (local_hour_3d.ge.(hh-1) .and. local_hour_3d.lt.hh)
      if (id_temp_hour(hh).gt.0) then
         used = send_data (id_temp_hour(hh),t,Time,is,js,1,mask=mask_local_hour)
                 used = send_data (id_temp_hour(hh),t,Time,is,js,1,mask=mask_local_hour)
      end if
      if (id_ps_hour(hh).gt.0) then
         used = send_data (id_ps_hour(hh),phalf(:,:,kd+1),Time,is_in=is,js_in=js,mask=mask_local_hour(:,:,1))
      end if
   end do

   do n=1,nt
      do hh=1,24
         mask_local_hour = (local_hour_3d.ge.(hh-1) .and. local_hour_3d.lt.hh)
         if (id_tracer_diag_hour(n,hh) .gt. 0 ) then
            used = send_data (id_tracer_diag_hour(n,hh),tracer_orig(:,:,:,n),Time,is,js,1,mask=mask_local_hour)
         end if
      end do

      if ( id_tracer_diag(n) .gt. 0 ) then

         call get_tracer_names (MODEL_ATMOS, n, name = tracer_name,  &
              units = tracer_units)
         if ( tracer_units .eq. "vmr" ) then
            used  = send_data (id_tracer_diag(n),     &
                 1.e3*rho(:,:,:)/WTMAIR * (tracer_orig(:,:,:,n)+rdt(:,:,:,n)), &
                 Time, is_in=is, js_in=js, ks_in=1)
         else
            used  = send_data (id_tracer_diag(n),     &
                 rho(:,:,:) * (tracer_orig(:,:,:,n)+rdt(:,:,:,n)), &
                 Time, is_in=is, js_in=js, ks_in=1)
         end if
      end if
   end do


 end subroutine atmos_tracer_driver
! </SUBROUTINE>

!#######################################################################

! <SUBROUTINE NAME="atmos_tracer_driver_init">
!   <OVERVIEW>
!     Subroutine to initialize the tracer driver module.
!   </OVERVIEW>
!   <DESCRIPTION>
!   The purpose of the arguments here are for passing on to the individual
!   tracer code. The user may wish to provide initial values which can be
!   implemented in the initialization part of the tracer code. Remember that
!   the tracer manager will provide a simple fixed or exponential profile if
!   the user provides data for this within the field table. However if a more
!   complicated profile is required then it should be set up in the
!   initialization section of the user tracer code.
!   </DESCRIPTION>
!   <TEMPLATE>
!     call atmos_tracer_driver_init (lonb,latb, r, mask, axes, Time)
!   </TEMPLATE>
!   <IN NAME="lonb" TYPE="real" DIM="(:,:)">
!     The longitude corners for the local domain.
!   </IN>
!   <IN NAME="latb" TYPE="real" DIM="(:,:)">
!     The latitude corners for the local domain.
!   </IN>
!   <IN NAME="mask" TYPE="real, optional" DIM="(:,:,:)">
!      optional mask (0. or 1.) that designates which grid points
!           are above (=1.) or below (=0.) the ground dimensioned as
!           (nlon,nlat,nlev).
!   </IN>
!   <IN NAME="Time" TYPE="type(time_type)">
!     Model time.
!   </IN>
!   <IN NAME="axes" TYPE="integer" DIM="(4)">
!     The axes relating to the tracer array dimensioned as
!      (nlon, nlat, nlev, ntime)
!   </IN>
!   <INOUT NAME="r" TYPE="real" DIM="(:,:,:,:)">
!     Tracer fields dimensioned as (nlon,nlat,nlev,ntrace).
!   </INOUT>
 subroutine atmos_tracer_driver_init (lonb, latb, r, axes, Time, phalf, mask)

!-----------------------------------------------------------------------
           real, intent(in),    dimension(:,:)             :: lonb, latb
           real, intent(inout), dimension(:,:,:,:)         :: r
type(time_type), intent(in)                                :: Time
        integer, intent(in)                                :: axes(4)
           real, intent(in),    dimension(:,:,:)           :: phalf
           real, intent(in),    dimension(:,:,:), optional :: mask

!-----------------------------------------------------------------------
! Local variables
!-----------------------------------------------------------------------
      integer :: nbr_layers
      integer :: unit, ierr, io, logunit, n, outunit
!<f1p
      character(len=32) :: tracer_units, tracer_name
      integer     :: hh
      character*4 :: hstr
      character(len=256) :: cmip_name,cmip_longname, cmip_longname2
      logical :: cmip_is_aerosol, do_pm, do_check
      real    :: tracer_mw, sum_N_ox
!>

!-----------------------------------------------------------------------
!
!  When initializing additional tracers, the user needs to make changes
!
!-----------------------------------------------------------------------

      if (module_is_initialized) return

!------------------------------------------------------------------------
!   read namelist.
!------------------------------------------------------------------------
      if ( file_exist('input.nml')) then
#ifdef INTERNAL_FILE_NML
        read (input_nml_file, nml=atmos_tracer_driver_nml, iostat=io)
        ierr = check_nml_error(io,'atmos_tracer_driver_nml')
#else
        unit =  open_namelist_file ( )
        ierr=1; do while (ierr /= 0)
        read  (unit, nml=atmos_tracer_driver_nml, iostat=io, end=10)
        ierr = check_nml_error(io, 'atmos_tracer_driver_nml')
        end do
10      call close_file (unit)
#endif
      endif
!--------- write version and namelist to standard log ------------
      call write_version_number (version, tagname)
      logunit=stdlog()
      if ( mpp_pe() == mpp_root_pe() ) &
                         write ( logunit, nml=atmos_tracer_driver_nml )

!---------------------------------------------------------------------
!  make sure that astronomy_mod has been initialized (if radiation
!  not being called in this run, it will not have previously been
!  initialized).
!---------------------------------------------------------------------
      call astronomy_init

!If we wish to automatically register diagnostics for wet and dry
! deposition, do it now.
      call atmos_tracer_utilities_init(lonb, latb, axes, Time)

!----- set initial value of radon ------------

      call atmos_radon_init(r, axes, Time, nradon, mask)
      radon_clock = mpp_clock_id( 'Tracer: Radon', &
           grain=CLOCK_MODULE )

!----- initialize the convection tracer ------------

      call atmos_convection_tracer_init(r, phalf, axes, Time, &
                                        nconvect,  mask)
      convect_clock = mpp_clock_id( 'Tracer: Convection tracer', &
           grain=CLOCK_MODULE )

!----- initialize the age tracer ------------

      call atmos_age_tracer_init( r, axes, Time, nage, lonb, latb, phalf, mask)
      age_tracer_clock = mpp_clock_id( 'Tracer: Age tracer', grain=CLOCK_MODULE )

      call get_number_tracers (MODEL_ATMOS, num_tracers=nt, &
                               num_prog=ntp)


!------------------------------------------------------------------------
! Initialize stratospheric chemistry
!------------------------------------------------------------------------
      nsphum   = get_tracer_index(MODEL_ATMOS,'sphum')
      no3      = get_tracer_index(MODEL_ATMOS,'O3')
      no3ch    = get_tracer_index(MODEL_ATMOS,'O3_chem')
      nextinct = get_tracer_index(MODEL_ATMOS,'Extinction')
      naerosol = get_tracer_index(MODEL_ATMOS,'Aerosol')
!
!  Set up photolysis rates etc.
!
      do_coupled_stratozone = strat_chem_driver_init()
      stratozone_clock = mpp_clock_id( 'Tracer: Stratospheric Ozone', &
           grain=CLOCK_MODULE )


!------------------------------------------------------------------------
! Initialize tropospheric chemistry
!------------------------------------------------------------------------
      allocate( drydep_data(nt) )
      do_tropchem = tropchem_driver_init(r,mask,axes,Time,lonb,latb,phalf,drydep_data)
      if ( .not. do_tropchem ) then
          do n = 1,ntp
             call dry_deposition_init(n,lonb,latb,drydep_data(n))
          end do
      end if
      tropchem_clock = mpp_clock_id( 'Tracer: Tropospheric chemistry', &
           grain=CLOCK_MODULE )

! ----------Interactive traceres--------------------
! If any of the interactive tracers are activated, get the
! tracer number and initialize it.
!      ntraceer= get_tracer_index(MODEL_ATMOS,'tracer_name')
!      if (ntracer > 0) then
!        call {tracer}_init(lonb, latb, r, axes, Time, mask)
!      endif


!    get tracer indices
      nbcphobic = get_tracer_index(MODEL_ATMOS,'bcphob')
      nbcphilic = get_tracer_index(MODEL_ATMOS,'bcphil')
      nomphobic = get_tracer_index(MODEL_ATMOS,'omphob')
      nomphilic = get_tracer_index(MODEL_ATMOS,'omphil')
      nsulfate=0
      nSO4      = get_tracer_index(MODEL_ATMOS,'simpleSO4')
      if (nSO4 .gt. 0) then
        nsulfate=nsulfate+1
        tr_nbr_sulfate(1)=nSO4
        do_tracer_sulfate(1)=.true.
      endif
      nSO2      = get_tracer_index(MODEL_ATMOS,'simpleSO2')
      if (nSO2 .gt. 0) then
        nsulfate=nsulfate+1
        tr_nbr_sulfate(2)=nSO2
        do_tracer_sulfate(2)=.true.
      endif
      nDMS      = get_tracer_index(MODEL_ATMOS,'simpleDMS')
      if (nDMS .gt. 0) then
        nsulfate=nsulfate+1
        tr_nbr_sulfate(3)=nDMS
        do_tracer_sulfate(3)=.true.
      endif
      nH2O2     = get_tracer_index(MODEL_ATMOS,'simpleH2O2')
      if (nH2O2 .gt. 0) then
        nsulfate=nsulfate+1
        tr_nbr_sulfate(4)=nH2O2
        do_tracer_sulfate(4)=.true.
      endif
      nMSA      = get_tracer_index(MODEL_ATMOS,'simpleMSA')
      if (nMSA .gt. 0) then
        nsulfate=nsulfate+1
        tr_nbr_sulfate(5)=nMSA
        do_tracer_sulfate(5)=.true.
      endif
      nNH4NO3   = get_tracer_index(MODEL_ATMOS,'nh4no3')
      nNH4      = get_tracer_index(MODEL_ATMOS,'nh4')
      nSOA      = get_tracer_index(MODEL_ATMOS,'SOA')
      nsf6      = get_tracer_index(MODEL_ATMOS,'sf6')
      nch3i     = get_tracer_index(MODEL_ATMOS,'ch3i')
      nco2      = get_tracer_index(MODEL_ATMOS,'co2')
      nDMS_cmip = get_tracer_index(MODEL_ATMOS,'DMS')
      nSO2_cmip = get_tracer_index(MODEL_ATMOS,'so2')
      nSO4_cmip = get_tracer_index(MODEL_ATMOS,'so4')
      nNH3_cmip = get_tracer_index(MODEL_ATMOS,'nh3')
      nOH       = get_tracer_index(MODEL_ATMOS,'oh')
      nC4H10    = get_tracer_index(MODEL_ATMOS,'c4h10')
! Check for presence of OH and C4H10 (diagnostic) tracers
! If not present set index to 1 so interface calls do not fail,
! but FATAL error will be issued by atmos_sulfate_init,
! atmos_carbon_aerosol_init (if do_dynamic_bc or do_dynamic_om), and
! atmos_SOA_init (if use_interactive_tracers)
      if (nOH == NO_TRACER) nOH = 1
      if (nC4H10 == NO_TRACER) nC4H10 = 1
      ncodirect = get_tracer_index(MODEL_ATMOS,'codirect')
      ne90      = get_tracer_index(MODEL_ATMOS,'e90')

! Number of vertical layers
      nbr_layers=size(r,3)

! initialize the tracers
!carbonaceous aerosols
      if (nbcphobic > 0 .or. nbcphilic >0 .or.  &
          nomphobic > 0 .or. nomphilic >0 ) then
        call atmos_carbon_aerosol_init(lonb, latb, axes, Time, mask)
        carbon_clock = mpp_clock_id( 'Tracer: Carbonaceous aerosol', &
                       grain=CLOCK_MODULE )

      endif
!dust aerosols
      call atmos_dust_init (lonb, latb, axes, Time, mask)
      if (do_dust) then
        dust_clock = mpp_clock_id( 'Tracer: Dust aerosol', &
                     grain=CLOCK_MODULE )
      endif
!sea salt
      call atmos_sea_salt_init (lonb, latb, axes, Time, mask)
      if (do_seasalt) then
        seasalt_clock = mpp_clock_id( 'Tracer: Seasalt aerosol', &
                     grain=CLOCK_MODULE )
      endif
!sulfur cycle
      if (nsulfate > 0 ) then
        call atmos_sulfate_init ( lonb, latb, nbr_layers, axes, Time, mask)
        sulfur_clock = mpp_clock_id( 'Tracer: Sulfur', &
                       grain=CLOCK_MODULE )
      endif
!SOA
      if ( nSOA > 0 ) then
        call atmos_SOA_init ( lonb, latb, nbr_layers, axes, Time, mask)
        SOA_clock = mpp_clock_id( 'Tracer: SOA', &
                    grain=CLOCK_MODULE )
      endif
!sf6
      if (nsf6 > 0) then
        call atmos_sulfur_hex_init (lonb, latb, r, axes, Time, mask)
        sf6_clock = mpp_clock_id( 'Tracer: SF6', &
                    grain=CLOCK_MODULE )
      endif
!ch3i
      if (nch3i > 0) then
        call atmos_ch3i_init (lonb, latb, axes, Time, mask)
        ch3i_clock = mpp_clock_id( 'Tracer: CH3I', &
                     grain=CLOCK_MODULE )
      endif

!co2
      if (nco2 > 0) then
      call atmos_co2_init ( Time, size(r,1), size(r,2), axes(1:3))
        co2_clock = mpp_clock_id( 'Tracer: CO2', &
                    grain=CLOCK_MODULE )
      endif

      if (ncodirect > 0 .or. ne90 > 0) then
        call regional_tracer_driver_init (lonb, latb, axes, Time, mask)
      endif

!tropopause diagnostics
      call atmos_tropopause_init (Time)
      tropopause_clock = mpp_clock_id( 'Tracer: Tropopause', grain=CLOCK_MODULE )

     call get_number_tracers (MODEL_ATMOS, num_tracers=nt, &
                               num_prog=ntp)

     id_landfr = register_diag_field ( mod_name,                    &
            'landfr_atm', axes(1:2), Time,               &
            'land fraction',                                 &
            'fraction', missing_value=-999.     )
     id_seaicefr = register_diag_field ( mod_name,                    &
            'seaicefr_atm', axes(1:2), Time,               &
            'seaice fraction',                                 &
            'fraction', missing_value=-999.     )
     id_snowfr = register_diag_field ( mod_name,                    &
            'snowfr_atm', axes(1:2), Time,               &
            'snow cover fraction',                                 &
            'fraction', missing_value=-999.     )
     id_vegnfr = register_diag_field ( mod_name,                    &
            'vegnfr_atm', axes(1:2), Time,               &
            'vegetation cover fraction',                                 &
            'fraction', missing_value=-999.     )
     id_vegnlai = register_diag_field ( mod_name,                    &
            'vegnlai_atm', axes(1:2), Time,               &
            'vegetation leaf area index',                                &
            'fraction', missing_value=-999.     )

      id_om_ddep = register_diag_field (mod_name, &
          'om_ddep', axes(1:2), Time, &
          'total dry deposition of om', 'kg/m2/s')

      id_bc_ddep = register_diag_field (mod_name, &
          'bc_ddep', axes(1:2), Time, &
          'total dry deposition of bc', 'kg/m2/s')

      id_nh4_ddep_cmip = register_diag_field (mod_name, &
          'tot_nh4_ddep_cmip', axes(1:2), Time, &
          'total dry deposition of ammonium', 'kg/m2/s')

      id_nh4no3_col = register_diag_field (mod_name, &
          'tot_no3_col', axes(1:2), Time, &
          'total aerosol load of nitrate', 'kg/m2')

      id_nh4_col  = register_diag_field (mod_name, &
          'tot_nh4_col', axes(1:2), Time, &
          'total aerosol load of ammonium', 'kg/m2')

      id_nh4no3_cmip = register_diag_field (mod_name, &
          'tot_no3', axes(1:3), Time, &
          'total nitrate', 'kg/m3')

      id_nh4_cmip  = register_diag_field (mod_name, &
          'tot_nh4', axes(1:3), Time, &
          'total ammonium', 'kg/m3')

      id_nh4no3_cmipv2 = register_diag_field (mod_name, &
          'tot_no3v2', axes(1:3), Time, &
          'total nitrate', 'kg/m3')

      id_nh4_cmipv2  = register_diag_field (mod_name, &
          'tot_nh4v2', axes(1:3), Time, &
          'total ammonium', 'kg/m3')

      id_so2_cmip  = register_diag_field (mod_name, &
          'so2_cmip', axes(1:3), Time, &
          'SO2', 'kg/m3')

      id_dms_cmip  = register_diag_field (mod_name, &
          'dms_cmip', axes(1:3), Time, &
          'DMS', 'kg/m3')

      id_so2_cmipv2  = register_diag_field (mod_name, &
          'so2_cmipv2', axes(1:3), Time, &
          'SO2', 'kg/m3')

      id_dms_cmipv2  = register_diag_field (mod_name, &
          'dms_cmipv2', axes(1:3), Time, &
          'DMS', 'kg/m3')

      !---- register cmip named variables ----

      ID_concso2 = register_cmip_diag_field_3d ( mod_name, &
                   'concso2', Time, 'Mole Fraction of SO2', '1.0', &
                   standard_name='mole_fraction_of_sulfur_dioxide_in_air')

      ID_concdms = register_cmip_diag_field_3d ( mod_name, &
                   'concdms', Time, 'Mole Fraction of DMS', '1.0', &
                    standard_name='mole_fraction_of_dimethyl_sulfide_in_air')

      if (nNH4NO3 > 0 .and. nNH4 > 0) then
        ID_concnh4 = register_cmip_diag_field_3d ( mod_name, &
                   'concnh4', Time, 'Concentration of NH4', 'kg m-3', &
                   standard_name='mass_concentration_of_ammonium_dry_aerosol_in_air')

        id_sconcnh4 = register_cmip_diag_field_2d ( mod_name, &
                  'sconcnh4', Time, 'Surface Concentration of NH4', 'kg m-3', &
                  standard_name='mass_concentration_of_ammonium_dry_aerosol_in_air')

        id_loadnh4 = register_cmip_diag_field_2d ( mod_name, &
                   'loadnh4', Time, 'Load of NH4', 'kg m-2', &
                   standard_name='atmosphere_mass_content_of_ammonium_dry_aerosol')

         id_drynh4 = register_cmip_diag_field_2d ( mod_name, &
                  'drynh4', Time, 'Dry Deposition Rate of NH4', 'kg m-2 s-1', &
                  standard_name='tendency_of_atmosphere_mass_content_of_ammonium_dry_aerosol_particles_due_to_dry_deposition')
      else
        id_sconcnh4 = 0
        id_loadnh4 = 0
        id_drynh4 = 0
      endif

      ID_airmass = register_cmip_diag_field_3d ( mod_name, &
           'airmass', Time, 'Vertically integrated mass content of air in layer', 'kg m-2', &
           standard_name='atmosphere_mass_of_air_per_unit_area')


      if (nNH4NO3 > 0) then
        ID_concno3 = register_cmip_diag_field_3d ( mod_name, &
                   'concno3', Time, 'Concentration of NO3 Aerosol', 'kg m-3', &
                   standard_name='mass_concentration_of_nitrate_dry_aerosol_in_air')

        id_sconcno3 = register_cmip_diag_field_2d ( mod_name, &
                  'sconcno3', Time, 'Surface Concentration of NO3', 'kg m-3', &
                  standard_name='mass_concentration_of_nitrate_dry_aerosol_in_air')

        id_loadno3 = register_cmip_diag_field_2d ( mod_name, &
                   'loadno3', Time, 'Load of NO3', 'kg m-2', &
                   standard_name='atmosphere_mass_content_of_nitrate_dry_aerosol')
      else
        id_sconcno3 = 0
        id_loadno3 = 0
      endif

      id_dryso2 = register_cmip_diag_field_2d ( mod_name, &
                     'dryso2', Time, 'Dry Deposition Rate of SO2', 'kg m-2 s-1', &
                     standard_name='tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_dry_deposition')

      id_dryso4 = register_cmip_diag_field_2d ( mod_name, &
                    'dryso4', Time, 'Dry Deposition Rate of SO4', 'kg m-2 s-1', &
                    standard_name='tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_dry_deposition')

      id_drydms = register_cmip_diag_field_2d ( mod_name, &
                    'drydms', Time, 'Dry Deposition Rate of DMS', 'kg m-2 s-1', &
                    standard_name='tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_dry_deposition')

      id_drynh3 = register_cmip_diag_field_2d ( mod_name, &
                    'drynh3', Time, 'Dry Deposition Rate of NH3', 'kg m-2 s-1', &
                    standard_name='tendency_of_atmosphere_mass_content_of_ammonia_due_to_dry_deposition')

      id_drybc = register_cmip_diag_field_2d ( mod_name, &
                  'drybc', Time, 'Dry Deposition Rate of Black Carbon Aerosol Mass', 'kg m-2 s-1', &
                  standard_name='tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_particles_due_to_dry_deposition')

      id_drypoa = register_cmip_diag_field_2d ( mod_name, &
                  'drypoa', Time, 'Dry Deposition Rate of Dry Aerosol Primary Organic Matter', 'kg m-2 s-1', &
                  standard_name='tendency_of_atmosphere_mass_content_of_primary_particulate_organic_matter_dry_aerosol_particles_due_to_dry_deposition')

      id_drysoa = register_cmip_diag_field_2d ( mod_name, &
                  'drysoa', Time, 'Dry Deposition Rate of Dry Aerosol Secondary Organic Matter', 'kg m-2 s-1', &
                  standard_name='tendency_of_atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol_particles_due_to_dry_deposition')

      id_dryoa = register_cmip_diag_field_2d ( mod_name, &
                  'dryoa', Time, 'Dry Deposition Rate of Dry Aerosol Total Organic Matter', 'kg m-2 s-1', &
                  standard_name='tendency_of_atmosphere_mass_content_of_particulate_organic_matter_dry_aerosol_particles_due_to_dry_deposition')
      !----

!<f1p: tracer diagnostics
      allocate( id_tracer_diag(nt) )
      allocate( id_tracer_diag_hour(nt,24) )
      id_tracer_diag(:) = 0

      allocate( ID_tracer_mol_mol(nt) )
      allocate( ID_tracer_kg_kg(nt) )
      allocate( id_tracer_col_kg_m2(nt) )
      allocate( id_tracer_ddep_kg_m2_s(nt) )
      allocate( id_tracer_surf_mol_mol(nt) )
      allocate( id_tracer_surf_kg_kg(nt) )

      allocate(conv_vmr_mmr(nt))
      conv_vmr_mmr = 1.
      allocate(nb_N(nt))
      allocate(nb_N_ox(nt))
      allocate(nb_N_red(nt))

      allocate(frac_pm1(nt))
      allocate(frac_pm10(nt))
      allocate(frac_pm25(nt))

      ID_BC = register_cmip_diag_field_3d ( mod_name, 'fam_bc_kg_kg', Time, &
             'Elemental carbon mass mixing ratio', 'kg kg-1', &
             standard_name='mass_fraction_of_elemental_carbon_dry_aerosol_particles_in_air')

      ID_OM = register_cmip_diag_field_3d ( mod_name, 'fam_om_kg_kg', Time, &
             'Total organic aerosol mass mixing ratio', 'kg kg-1', &
             standard_name='mass_fraction_of_particulate_organic_matter_dry_aerosol_particles_in_air')

      ID_DUST = register_cmip_diag_field_3d ( mod_name, 'fam_dust_kg_kg', Time, &
             'Dust aerosol mass mixing ratio', 'kg kg-1', &
             standard_name='mass_fraction_of_dust_dry_aerosol_particles_in_air')

      ID_SS = register_cmip_diag_field_3d ( mod_name, 'fam_seasalt_kg_kg', Time, &
             'Sea salt mass mixing ratio', 'kg kg-1', &
             standard_name='mass_fraction_of_seasalt_dry_aerosol_particles_in_air')

      ID_pm10 = register_cmip_diag_field_3d ( mod_name, 'pm10_kg_kg', Time, &
              'PM10 mass mixing ratio',  'kg kg-1', &
              standard_name='mass_fraction_of_pm10_ambient_aerosol_particles_in_air')

      ID_pm25 = register_cmip_diag_field_3d ( mod_name, 'pm25_kg_kg', Time, &
              'PM2.5 mass mixing ratio', 'kg kg-1', &
              standard_name='mass_fraction_of_pm2p5_dry_aerosol_particles_in_air')

      ID_pm1  = register_cmip_diag_field_3d ( mod_name, 'pm1_kg_kg', Time, &
              'PM1.0 mass mixing ratio', 'kg kg-1', &
              standard_name='mass_fraction_of_pm1_dry_aerosol_particles_in_air')

      id_pm25_surf = register_cmip_diag_field_2d ( mod_name, 'pm25_surf_kg_kg', Time, &
              'PM2.5 mass mixing ratio in lowest model layer', 'kg kg-1', &
              standard_name='mass_fraction_of_pm2p5_ambient_aerosol_particles_in_air')

      do_pm = .false.
      if (query_cmip_diag_id(ID_pm10) .or. query_cmip_diag_id(ID_pm1) .or. &
          query_cmip_diag_id(ID_pm25) .or. id_pm25_surf > 0) do_pm = .true.

      id_bc_col_kg_m2 = register_cmip_diag_field_2d ( mod_name, 'fam_bc_col_kg_m2', &
                       Time, 'Load of Black Carbon Aerosol', 'kg m-2', &
                       standard_name='atmosphere_mass_content_of_black_carbon_dry_aerosol')

      id_dust_col_kg_m2 = register_cmip_diag_field_2d ( mod_name, 'fam_dust_col_kg_m2', &
                         Time, 'Load of Dust', 'kg m-2', &
                         standard_name='atmosphere_mass_content_of_dust_dry_aerosol')

      id_seasalt_col_kg_m2 = register_cmip_diag_field_2d ( mod_name, 'fam_seasalt_col_kg_m2', &
                            Time, 'Load of Seasalt', 'kg m-2', &
                            standard_name='atmosphere_mass_content_of_seasalt_dry_aerosol')

      id_toz = register_cmip_diag_field_2d ( mod_name, 'toz', &
                       Time, 'Total Ozone Column', 'm', &
                       standard_name='equivalent_thickness_at_stp_of_atmosphere_ozone_content')

      id_tropoz = register_cmip_diag_field_2d ( mod_name, 'tropoz', &
                       Time, 'Tropospheric Ozone Column', 'm', &
                       standard_name='equivalent_thickness_at_stp_of_atmosphere_ozone_content')

      outunit = stdout()

      do n = 1,nt
         call get_tracer_names (MODEL_ATMOS, n, name = tracer_name,  &
              units = tracer_units)

         if ( tracer_units .eq. "vmr" ) then
            id_tracer_diag(n)  = register_diag_field (mod_name, &
                 tracer_name, axes(1:3), Time, &
                 tracer_name, 'mole/m3')
         elseif ( tracer_units .eq. "mmr" ) then
            id_tracer_diag(n)  = register_diag_field (mod_name, &
                 tracer_name, axes(1:3), Time, &
                 tracer_name, 'kg/m3')
         end if

         do hh=1,24
            write(hstr,'(A1,I2.2,A1)') '_',hh-1,'h'
            id_tracer_diag_hour(n,hh)  = register_diag_field (mod_name, &
                 trim(tracer_name)//hstr, axes(1:3), Time, &
                 trim(tracer_name)//hstr, tracer_units, missing_value=-999.,mask_variant = .true.)
         end do

         call  get_cmip_param (n, cmip_name=cmip_name, cmip_longname=cmip_longname, cmip_longname2=cmip_longname2)
         call  get_chem_param (n, mw=tracer_mw, conv_vmr_mmr=conv_vmr_mmr(n), is_aerosol=cmip_is_aerosol, &
                               nb_N=nb_N(n), nb_N_Ox=nb_N_Ox(n), nb_N_red=nb_N_red(n), &
                               frac_pm1=frac_pm1(n), frac_pm25=frac_pm25(n), frac_pm10=frac_pm10(n))

         write(outunit,'(a,i3)') 'n=',n
         write(outunit,'(7a)') 'tracer_name="',trim(tracer_name),'", cmip_name="',trim(cmip_name),'", cmip_longname="',trim(cmip_longname),'"'
         write(outunit,'(5(a,g14.6))') 'mwt=',tracer_mw, ', conv_vmr_mmr=',conv_vmr_mmr(n), &
                                       ', nb_N=',nb_N(n),', nb_N_ox=',nb_N_ox(n),', nb_N_red=',nb_N_red(n)
         write(outunit,'(3(a,f7.4))') 'frac_pm1=',frac_pm1(n), ', frac_pm25=',frac_pm25(n), ', frac_pm10=',frac_pm10(n)

         ID_tracer_mol_mol(n) = register_cmip_diag_field_3d ( mod_name, &
              trim(tracer_name)//'_mol_mol', Time, &
              trim(cmip_longname)//' volume mixing ratio', 'mol mol-1', &
              standard_name='mole_fraction_of_'//trim(cmip_name)//'_in_air')

         ID_tracer_kg_kg(n)   = register_cmip_diag_field_3d ( mod_name, &
              trim(tracer_name)//'_kg_kg', Time, &
              'Aerosol '//trim(cmip_longname)//' mass mixing ratio', 'kg kg-1', &
              standard_name='mass_fraction_of_'//trim(cmip_name)//'_dry_aerosol_particles_in_air')

         id_tracer_surf_mol_mol(n) = register_cmip_diag_field_2d ( mod_name, &
              trim(tracer_name)//'_surf_mol_mol', Time, &
              trim(cmip_longname)//' volume mixing ratio in lowest model layer', 'mol mol-1', &
              standard_name='mole_fraction_of_'//trim(cmip_name)//'_in_air')

         id_tracer_col_kg_m2(n) = register_cmip_diag_field_2d ( mod_name, &
              trim(tracer_name)//'_col_kg_m2', Time, &
              "Load of "//trim(cmip_longname2), 'kg m-2', &
              standard_name='atmosphere_mass_content_of_'//trim(cmip_name)//'_dry_aerosol')


         id_tracer_surf_kg_kg(n) = register_cmip_diag_field_2d ( mod_name, &
              trim(tracer_name)//'_surf_kg_kg', Time, &
              trim(cmip_longname)//' mass mixing ratio in lowest model layer', 'kg kg-1', &
              standard_name='mass_fraction_of_'//trim(cmip_name)//'_in_air')


         if (cmip_is_aerosol) then
            id_tracer_ddep_kg_m2_s(n) = register_cmip_diag_field_2d ( mod_name, &
                 trim(tracer_name)//'_ddep_kg_m2_s', Time, 'Dry Deposition Rate of '//trim(cmip_longname), 'kg m-2 s-1', &
                 standard_name='tendency_of_atmosphere_mass_content_of_'//trim(cmip_name)//'_dry_aerosol_particles_due_to_dry_deposition')
         else
            id_tracer_ddep_kg_m2_s(n) = register_cmip_diag_field_2d ( mod_name, &
                 trim(tracer_name)//'_ddep_kg_m2_s', Time, 'Dry Deposition Rate of '//trim(cmip_longname), 'kg m-2 s-1', &
                 standard_name='tendency_of_atmosphere_mass_content_of_'//trim(cmip_name)//'_due_to_dry_deposition')
         end if

         ! sanity check
         do_check = .false.
         if (id_tracer_ddep_kg_m2_s(n) > 0 .or. id_tracer_surf_kg_kg(n) > 0 .or. &
             id_tracer_col_kg_m2(n) > 0 .or.  query_cmip_diag_id(ID_tracer_kg_kg(n))) do_check = .true.
         if (do_pm .and. (frac_pm1(n) > 0.0 .or. frac_pm10(n) > 0.0 .or. frac_pm25(n) > 0.0)) do_check = .true.
         if (do_check .and. conv_vmr_mmr(n) < 0.0) then
            call error_mesg ('Tracer_driver', 'mw needs to be defined for tracer: '//trim(tracer_name), FATAL)
         end if
      end do

      do hh=1,24
         write(hstr,'(A1,I2.2,A1)') '_',hh-1,'h'
         id_temp_hour(hh)  = register_diag_field (mod_name, &
              'temp'//hstr, axes(1:3), Time, &
              'temp'//hstr, 'K', missing_value=-999.,mask_variant = .true.)
         id_ps_hour(hh)  = register_diag_field (mod_name, &
              'ps'//hstr, axes(1:2), Time, &
              'ps'//hstr, 'Pa', missing_value=-999.,mask_variant = .true.)
      end do

      id_n_ox_ddep = 0
      id_n_ddep = 0
      id_n_red_ddep=0

      write (outunit,*) 'fam_N_ox is comprised of :'
      sum_N_ox = 0.0
      do n = 1,nt
         if ( nb_N_ox(n) .gt. 0.) then
            sum_N_ox = sum_N_ox + nb_N_ox(n)
            call get_tracer_names (MODEL_ATMOS, n, name = tracer_name, units = tracer_units)
            write (outunit,'(2a,g14.6)') trim(tracer_name),', nb_N_ox=',nb_N_ox(n)
         end if
      end do
      if (sum_N_ox > 0.0) then
        id_n_ox_ddep =  register_cmip_diag_field_2d ( mod_name, &
                'fam_noy_ddep_kg_m2_s', Time, &
                'Dry Deposition Rate of all Nitrogen Oxides (NOY)', 'kg m-2 s-1', &
                standard_name='tendency_of_atmosphere_mass_content_of_noy_expressed_as_nitrogen_due_to_dry_deposition')
      end if

      write (outunit,*) 'fam_N_red is comprised of :'
      do n = 1,nt
         if ( nb_N_red(n) .gt. 0.) then
            call get_tracer_names (MODEL_ATMOS, n, name = tracer_name, units = tracer_units)
            write (outunit,'(2a,g14.6)') trim(tracer_name),', nb_N_red=',nb_N_red(n)
         end if
      end do

      write (outunit,*) 'fam_N is comprised of :'
      do n = 1,nt
         if ( nb_N(n) .gt. 0.) then
            call get_tracer_names (MODEL_ATMOS, n, name = tracer_name, units = tracer_units)
            write (outunit,'(2a,g14.6)') trim(tracer_name),', nb_N=',nb_N(n)
         end if
      end do

  !BW write (outunit,*) 'frac_pm25 is comprised of :'
      write (outunit,*) 'pm** is comprised of :'
      do n = 1,nt
         if ( nb_N(n) .gt. 0.) then
            call get_tracer_names (MODEL_ATMOS, n, name = tracer_name, units = tracer_units)
  !BW       write (outunit,'(2a,f7.4)') trim(tracer_name),', frac_pm25=',frac_pm25(n)
            write (outunit,'(a,3(a,f7.4))') trim(tracer_name),', frac_pm25=',frac_pm25(n), &
                                            ', frac_pm10=',frac_pm10(n),', frac_pm1=',frac_pm1(n)
         end if
      end do

  !BW write (outunit,*) 'frac_pm10 is comprised of :'
  !BW do n = 1,nt
  !BW    if ( nb_N(n) .gt. 0.) then
  !BW       call get_tracer_names (MODEL_ATMOS, n, name = tracer_name, units = tracer_units)
  !BW       write (outunit,'(2a,f7.4)') trim(tracer_name),' frac_pm10=',frac_pm10(n)
  !BW    end if
  !BW end do

  !BW write (outunit,*) 'frac_pm1 is comprised of :'
  !BW do n = 1,nt
  !BW    if ( nb_N(n) .gt. 0.) then
  !BW       call get_tracer_names (MODEL_ATMOS, n, name = tracer_name, units = tracer_units)
  !BW       write (outunit,'(2a,f7.4)') trim(tracer_name),' frac_pm1=',frac_pm1(n)
  !BW    end if
  !BW end do
!>

      module_is_initialized = .TRUE.

 end subroutine atmos_tracer_driver_init



!#####################################################################

subroutine atmos_tracer_driver_time_vary (Time)

type(time_type), intent(in) :: Time

      if (nbcphobic > 0 .and. nbcphilic > 0 .and. &
          nomphobic > 0 .and. nomphilic > 0) then
        call atmos_carbon_aerosol_time_vary (Time)
      endif

      if (nch3i > 0) then
        call atmos_ch3i_time_vary (Time)
      endif

      if (do_dust) then
        call atmos_dust_time_vary (Time)
      endif

      if (nSOA > 0 ) then
        call atmos_SOA_time_vary (Time)
      endif

      if (nDMS > 0 .or. nSO2 > 0 .or. nSO4 > 0 &
                   .or. nMSA > 0 .or. nH2O2 > 0 ) then
        call atmos_sulfate_time_vary (Time)
      endif

      if (nco2 > 0) then
        call atmos_co2_time_vary (Time)
      endif

      call dry_deposition_time_vary (drydep_data, Time)

      if (do_tropchem) then
        call tropchem_driver_time_vary (Time)
      endif



end subroutine atmos_tracer_driver_time_vary


!#####################################################################

subroutine atmos_tracer_driver_endts

      if (do_tropchem) then
        call tropchem_driver_endts
      endif
      if (nbcphobic > 0 .and. nbcphilic > 0 .and. &
          nomphobic > 0 .and. nomphilic > 0) then
        call atmos_carbon_aerosol_endts
      endif
      if (nch3i > 0) then
        call atmos_ch3i_endts
      endif
      if (do_dust) then
        call atmos_dust_endts
      endif
      if (nSOA > 0 ) then
        call atmos_soa_endts
      endif
      if (nDMS > 0 .or. nSO2 > 0 .or. nSO4 > 0 &
                   .or. nMSA > 0 .or. nH2O2 > 0 ) then
        call atmos_sulfate_endts
      endif
      call dry_deposition_endts (drydep_data)

end subroutine atmos_tracer_driver_endts




!#####################################################################


! </SUBROUTINE>

!#######################################################################

! <SUBROUTINE NAME="atmos_tracer_driver_end">
!   <OVERVIEW>
!     Subroutine to terminate the tracer driver module.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Termination routine for tracer_driver. It should also call
!     the destructors for the individual tracer routines.
!   </DESCRIPTION>
!   <TEMPLATE>
!     call atmos_tracer_driver_end
!   </TEMPLATE>
 subroutine atmos_tracer_driver_end

!-----------------------------------------------------------------------
integer :: logunit
      if (mpp_pe() /= mpp_root_pe()) return

      logunit=stdlog()
      write (logunit,'(/,(a))') 'Exiting tracer_driver, have a nice day ...'

      call atmos_radon_end
      call atmos_sulfur_hex_end
      call atmos_convection_tracer_end
      call atmos_dust_end
      call atmos_sea_salt_end
      call atmos_sulfate_end
      call atmos_SOA_end
      call atmos_carbon_aerosol_end
      if (nch3i > 0) then
        call atmos_ch3i_end
      endif
      call atmos_age_tracer_end
      call atmos_co2_end

!for cmip6 (f1p)
      deallocate( id_tracer_diag )
      deallocate( ID_tracer_mol_mol )
      deallocate( ID_tracer_kg_kg )
      deallocate( id_tracer_surf_mol_mol )
      deallocate( id_tracer_col_kg_m2 )
      deallocate( id_tracer_surf_kg_kg )
      deallocate( conv_vmr_mmr )
      deallocate( id_tracer_ddep_kg_m2_s )

      deallocate(nb_N_red)
      deallocate(nb_N_ox)
      deallocate(nb_N)
      deallocate(frac_pm1)
      deallocate(frac_pm25)
      deallocate(frac_pm10)
      deallocate( id_tracer_diag_hour )

      module_is_initialized = .FALSE.

!-----------------------------------------------------------------------

 end subroutine atmos_tracer_driver_end
! </SUBROUTINE>

!#######################################################################
! <SUBROUTINE NAME="atmos_tracer_flux_init">
!   <OVERVIEW>
!     Subroutine to initialize the ocean-atmosphere gas flux modules
!   </OVERVIEW>
!   <DESCRIPTION>
!     Subroutine to initialize the ocean-atmosphere gas flux modules
!   </DESCRIPTION>

subroutine atmos_tracer_flux_init

call atmos_co2_flux_init
call atmos_dust_flux_init

return

end subroutine atmos_tracer_flux_init
! </SUBROUTINE>

!######################################################################
! <SUBROUTINE NAME="atmos_tracer_driver_gather_data">
!   <OVERVIEW>
!     Subroutine to terminate the tracer driver module.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Termination routine for tracer_driver. It should also call
!     the destructors for the individual tracer routines.
!   </DESCRIPTION>
!   <TEMPLATE>
!     call atmos_tracer_driver_gather_data
!   </TEMPLATE>
 subroutine atmos_tracer_driver_gather_data(gas_fields, tr_bot)

use coupler_types_mod, only: coupler_2d_bc_type

type(coupler_2d_bc_type), intent(inout) :: gas_fields
real, dimension(:,:,:), intent(in)      :: tr_bot

!-----------------------------------------------------------------------

  call atmos_co2_gather_data(gas_fields, tr_bot)

!-----------------------------------------------------------------------

 end subroutine atmos_tracer_driver_gather_data

 subroutine atmos_tracer_driver_gather_data_down(gas_fields, tr_bot)

use coupler_types_mod, only: coupler_2d_bc_type

type(coupler_2d_bc_type), intent(inout) :: gas_fields
real, dimension(:,:,:), intent(in)      :: tr_bot

!-----------------------------------------------------------------------

  call atmos_dust_gather_data(gas_fields, tr_bot)

!-----------------------------------------------------------------------

 end subroutine atmos_tracer_driver_gather_data_down
! </SUBROUTINE>
!######################################################################
! given a tracer index, returns true if this tracer has non-zero
! sedimentation flux at the bottom of the atmosphere
function atmos_tracer_has_surf_setl_flux(tr) result(ret)
   logical :: ret
   integer, intent(in) :: tr ! tracer index

   ret=.FALSE.
   if(is_dust_tracer(tr)) then
       ret = dust_has_surf_setl_flux(tr)
   ! later we can add:
   ! elseif (is_seasalt_tracer(tr)) then
   !   ret = salt_has_surf_setl_flux(tr)
   ! etc...
   endif
end function atmos_tracer_has_surf_setl_flux

!######################################################################
! given a tracer index, returns sedimentation flux at the bottom of
! the atmosphere for this tracer
subroutine get_atmos_tracer_surf_setl_flux(tr, setl_flux, dsetl_dtr)
  integer, intent(in)  :: tr         ! tracer index
  real,    intent(out) :: setl_flux(:,:) ! sedimentation flux at the bottom of the atmosphere
  real,    intent(out) :: dsetl_dtr(:,:) ! derivative of sedimentation flux w.r.t.
       ! the tracer concentration in the bottom layer

  setl_flux(:,:) = 0.0 ; dsetl_dtr(:,:) = 0.0
  if (is_dust_tracer(tr)) &
     call get_dust_surf_setl_flux(tr, setl_flux, dsetl_dtr)
  !if (is_seasalt_tracer(tr)) &
  !   call get_seasalt_surf_setl_flux(tr, setl_flux, dsetl_dtr)
end subroutine get_atmos_tracer_surf_setl_flux
!######################################################################


end module atmos_tracer_driver_mod
